Model construction and superconductivity analysis of organic conductors \beta-(BDA-TTP)_2MF_6 (M=P, As, Sb, Ta) based on first principles band calculation

We perform a first principles band calculation for a group of quasi-two-dimensional organic conductors \beta-(BDA-TTP)2MF6 (M=P, As, Sb, Ta). The ab-initio calculation shows that the density of states (DOS) is correlated with the band width of singly occupied (highest) molecular orbital (SOMO), while it is not necessarily correlated with the unit cell volume. The direction of the major axis of the cross section of the Fermi surface lies in the \Gamma-B direction, which differs from that obtained by the extended Huckel calculation. Then, we construct a tight-binding model which accurately reproduces the ab-initio band structure. The obtained transfer energies give smaller dimerization than in the extended Huckel band. As for the difference of the anisotropy of the Fermi surface, the transfer energies along the inter-stacking direction are smaller than those obtained in the extended Huckel calculation. Assuming spin-fluctuation-mediated superconductivity, we apply random phase approximation (RPA) to a two-band Hubbard model. This two-band Hubbard model is composed of the tight-binding model derived from the first principles band structure and an on-site (intra-molecule) repulsive interaction taken as a variable parameter. The obtained superconducting gap changes sign four times along the Fermi surface like in a d-wave gap, and the nodal direction is different from that obtained in the extended Huckel model. Anion dependence of Tc is qualitatively consistent with the experimental observation.


Introduction
Ever since their discovery, physics and chemistry of organic superconductors have attracted much attention [1]. One of the interesting features of these superconductors is the possibility of unconventional pairing arising from electron correlation effects. From this point of view, there have been efforts to control the strength of the electron correlation via chemical modification of molecules, and this has lead to the discovery of several new superconductors [2]. Among those superconductors, here we focus on a group of organic conductors β-(BDA-TTP) 2 MF 6 (M =P, As, Sb, Ta), where BDA-TTP is an abbreviation of 2,5-bis(1,3-dithian-2-ylidene)-1,3,4,6tetrathiapentalene. BDA-TTP molecule is derived from BDH-TTP by stereochemical modification, where BDH-TTP is an abbreviation of 2,5-bis(1,3-dithiolan-2-ylidene)-1,3,4,6-tetrathiapentalene. While BDH-TTP salts show stable metallic behavior suggesting weak electron correlation, the BDA-TTP salts exhibit stronger electron correlation, which is expected from the band width narrowing due to the chemical modification. In fact, β-(BDA-TTP) 2 MF 6 exhibit superconductivity with T c =7. 5, 5.8, 5.9K for M =Sb, As, P, respectively [3], while M =Ta does not superconduct [4]. The conducting layer is sandwiched by insulating layers of anion MF 6 as shown in figure  1(a), from which we may expect these salts to be quasi-two-dimensional. The alignment of the BDA-TTP donors in the conducting layer takes the β-type configuration shown in figure 1(b). Measurements of Shubnikov-de Haas (SdH) oscillation and angular- dependent magneto-resistance oscillation (AMRO) for β-(BDA-TTP) 2 SbF 6 have shown a two-dimensional Fermi surface [5] , which is consistent with the tight-binding band structure [3] obtained from the extended Hückel calculation [6]. However, a very recent AMRO measurement shows that the direction of the anisotropy of the Fermi surface differs from the previous AMRO measurement, and suggests the direction of the major axis of the cross section lies in the Γ-B direction of the Brillouin zone (see figure 2) [7]. It is worth mentioning that the effective cyclotron mass takes a large value of m * c /m 0 = 12.4±1.1, which indicates that these materials are strongly correlated electron systems, although the normal state exhibits a metallic nature.
Measurements of the upper critical field (H c2 ) of β-(BDA-TTP) 2 SbF 6 have shown spin-singlet pairing [8,9]. The temperature dependence of the electronic specific heat has indicated presence of a nodal structure in the superconducting order parameter [8], and a scanning tunneling microscope (STM) experiment also suggests presence of anisotropy in the gap, where the nodes (or the gap minima) lie in the k a ± k c direction [10]. On the other hand, a recent study of the in-plane anisotropy of H c2 shows that the H c2 maxima appears when the magnetic field is applied in the k a ± k c direction, i.e., the nodal direction indicated in the STM measurement [11]. This is puzzling because the H c2 maxima usually appears when the field is applied in the antinodal direction. A uniaxial pressure experiment for M =As and Sb shows that T c once increases and takes a maximum under the compression parallel to the stacking direction [12].
Theoretically, several studies have adopted the models derived from extended Hückel calculation, and analysis on superconductivity has also been performed assuming spin fluctuation mediated pairing [13,14]. Nonoyama et al. [15,16] applied RPA to the two band model for M =As and Sb, while Suzuki et al. [17] applied the fluctuation exchange (FLEX) approximation to the original two-band model and the single-band dimer model for M =As and Sb, and discussed the pressure dependence. There, the d-wave like gap function has been obtained, whose nodal direction is consistent with the STM measurement [10], but not with the recent study of in-plane anisotropy of H c2 [11].
Thus, there are several controversies regarding the relation between experiments and theoretical studies based on the extended Hückel model. In the present paper, we study theoretically the organic conductors β-(BDA-TTP) 2 MF 6 (M =P, As, Sb, Ta) based on a model constructed from ab-initio band calculation based on the density functional theory (DFT) exploiting the WIEN2K package [18]. In fact, the importance of model construction of molecular solids based on ab-initio calculation has recently pointed out in several studies [19,20,21]. We point out the difference in the band structure and the Fermi surface between the ab-initio (non-empirical method) and the extended Hückel calculations (semi-empirical method), and its consequence to the superconducting gap mediated by spin fluctuations. We also study the anion dependence of the band structure, and how it can affect the superconducting transition temperature.

ab-initio band calculation and model construction
We perform an ab-initio band calculation using all-electron full potential linearized augmented plane-wave (LAPW) + local orbitals (lo) method within WIEN2K [18]. This implements the DFT with different possible approximation for the exchange correlation potentials. The exchange correlation potential is calculated using the generalized gradient approximation (GGA).
To attain convergence in the eigenvalue calculation, the single-particle wave functions in the interstitial region are expanded by plane waves with a cut-off of R MT K max = 3.0, where R MT denotes the smallest muffin tin radius and K max is the maximum value of K vector in the plane wave expansion. For β-(BDA-TTP) 2 SbF 6 , for instance, the muffin-tin radii are taken as 1.74, 1.74, 1.62, 0.83, and 0.45 in atomic units (au) for Sb, F, S, C, and H, respectively. Thus K max is 3/0.45≈6.7, and the plane wave cutoff energy is 604.7 eV. Calculations are performed using 7×3×9 k-points in the irreducible Brillouin zone. We adopt the lattice structure determined experimentally (including the coordinates of the H atoms) [3], and we do not relax the atomic positions in the calculation.
Having done the ab-initio calculation, we then construct a tight-binding model which accurately reproduces the ab-initio band structure. From figure 1(b), it can be seen that the conducting BDA-TTP layer has two BDA-TTP molecules in a unit cell, so we regard one molecule as a site, and consider a two-band (two sites per unit cell) tight-binding model to fit the ab-initio band structure. The tight-binding Hamiltonian we consider, H 0 , is written in the form where i and j are unit cell indices, α and β specifies the sites in a unit cell, c † iασ (c iασ ) is a creation (annihilation) operator with spin σ at site α in the i-th unit cell, t iα:jβ is the electron transfer energy between (i, α) site and (j, β) site, and iα : jβ represents the summation over the bonds corresponding to the transfer.
By Fourier transformation, equation (1) is rewritten as where ε αβ (k ) is the site-indexed kinetic energy represented in k -space. The band dispersion is given by diagonalizing ε αβ (k ), where ξ γ (k ) is the dispersion of the γ band measured from the chemical potential, and d αγ (k ) is the unitary matrix that gives the transformation. In order to study unconventional superconductivity, we employ the two-band Hubbard model obtained by adding the on-site (intra-molecule) repulsive interaction to the tight-binding model derived from the ab-initio band structure. The Hubbard Hamiltonian, H, is where U is the on-site interaction and n iασ is the number operator of the electron on α-site in the i-th unit cell. The on-site U in the present study is fixed at a certain value. Although this choice does not have a quantitative basis, taking other values of U does not affect qualitatively the conclusion of the present paper regarding the superconducting gap symmetry and the anion dependence of T c . Note also that we ignore the off-site (inter-molecule) repulsive interactions. This is based on the result obtained in a previous theoretical work for the extended Hückel model of the present materials [15]. There it has been shown that the introduction of moderate inter-molecule interactions only slightly modifies the position of the nodes of the gap. Although the model adopted here is different, we believe that the effect of the inter-site interaction should be similar as far as the off-site interaction values are moderate.

Random phase approximation and superconducting gap equation
We apply RPA to the two-band Hubbard model given by equation (4) as follows. The bare susceptibility matrix in the site-representation is given by where N is the total number of unit cells, and f (ξ) is the Fermi distribution function (this is the multi-site version of the Lindhard function). The indices αβ on the left hand side means (α, β)-element of the bare susceptibility matrixX 0 . Within RPA, the spin and charge susceptibilities are obtained aŝ respectively.X sp ,X ch ,X 0 are the spin, charge, and bare susceptibility matrices, respectively,Û is the on-site interaction matrix, andÎ is the unit matrix. In a twoband system,Û ,X 0 ,X sp andX ch all become 2×2 matrices. The "spin susceptibility" (see figure 7(a)) in the present paper is obtained as the larger eigenvalue of the 2×2 spin susceptibility matrix. By using these susceptibility matrices, the pairing interactionV singlet for the spinsinglet state is given aŝ Using the pairing interaction given by equation (8), we solve the linearized gap equation to obtain the transition temperature T c and the superconducting gap function for the spin-singlet pairing state. The linearized gap equation within the weak-coupling theory is given by where λ is the eigenvalue of the linearized gap equation, and ϕ αβ (k ) is the (α, β)element of the gap function matrix. The transition temperature T c is the temperature where λ reaches unity. In a two-band system,V singlet andφ become 2×2 matrices. Although RPA is quantitatively insufficient for evaluating the absolute value of T c , we expect this approach to be valid for (i) studying the form of the superconducting gap function assuming spin-fluctuation-mediated pairing, and (ii) qualitatively comparing the tendency toward superconductivity (or T c ) between the different anions. We take 128×128 k-point meshes in the RPA calculation throughout the entire study.

ab-initio band calculation
We show in figure 2(a) the ab-initio band structure for β-(BDA-TTP) 2 SbF 6 at ambient pressure and room temperature. It can be seen that the singly occupied (highest) molecular orbital (SOMO) and the highest-occupied molecular orbital (HOMO) near the Fermi level (zero energy) are isolated from other bands, although there are two bands somewhat close to the HOMO band, which correspond to HOMO−1 and HOMO−2. Density of states (DOS) in figure 2(b) also shows the isolation of the HOMO and SOMO from the other bands. From these results, it can be considered that the SOMO and HOMO bands are the target bands to construct a low energy tight-binding model. figure 2(c) shows the Fermi surface of the ab-initio calculation, where the name of the k-points are presented only on the k Y (k b ) = 0 plane, e.g Γ=(0, 0, 0), Z=(0, 0, π/c ) and B=(π/a, 0, −π/c ). We also show similar plots for β-(BDA-TTP) 2 AsF 6 at ambient pressure and room temperature in figure 3. The band structure and the Fermi surface are similar between the two salts. M =P and Ta also give similar results (not shown). The DOS also looks similar qualitatively, but its quantitative dependence on the anions will be discussed later.
The Fermi surface is clearly cylindrical, reflecting the strong two dimensionality of this salt. The direction of the anisotropy of the Fermi surface differs from that found in the previous extended Hückel calculation [3,4,15,17]. The origin of this difference is intimately related to including/excluding the 3d-orbital of the sulfur atoms of the TTP skeleton donors in the extended Hückel calculation [22,23,24]. Actually, the shape of the Fermi surface given by the extended Hückel calculation in which the 3d-orbital of the sulfur atoms is excluded [25] is similar to the present ab-initio result. Focusing on the comparison of the experiments with our results of the Fermi surface, the anisotropy of the Fermi surface differs from the previous experiment [5], but a very recent AMRO measurement shows that the long-axis of the Fermi surface ellipse is in the Γ-B direction, which is at least qualitatively consistent with the present result [7]. However, the size of the Fermi surface is not consistent between this experiment and the present calculation, and this remains as a puzzle. Now, one of the important issues we study in the present paper is the effect of substituting the anion. Here we consider the anion dependence of the band-width and the DOS at the Fermi level for β-(BDA-TTP) 2 MF 6 (M=P, As, Sb, Ta) at ambient pressure and room temperature, and also for β-(BDA-TTP) 2 SbF 6 at 12K. As for the band width, we introduce W SH and W SOMO , where the former is the energy difference between the top of SOMO and the bottom of HOMO, and the latter is the band width of SOMO alone. figure 4(a) shows that W SH becomes smaller as the anion is enlarged from M =P to Sb, resulting in an increase of the unit-cell volume. The trend can naturally be understood as the effect of the chemical pressure. M =Ta, however, is an exception, where the increase of the unit-cell volume from M =Sb results in an increase in the band width. The relation between M =Sb and M =Sb(LT) shows that lowering the temperature results in an increases of W SH due to thermal compression. We show in (a) SOMO+HOMO band width, (b) SOMO band width, and (c) the density of states against the unit-cell volume. LT in parentheses for Sb means that the low temperature lattice parameters are used. The inset of (b) shows the schematic relation between SOMO and HOMO. figure 4(b) W SOMO against the unit-cell volume. The correlation between W SOMO and the volume is not monotonic in contrast to the volume dependence of W SH , even within the pnictogens M =P, As and Sb. The reason for this lies in the small dimerization gap between the SOMO and HOMO bands. Even though W SOMO of M =As is larger than that of M =P, W SH of As is smaller because the dimerization gap of As is smaller than that of P. A schematic figure is shown in the inset of figure 4(b), where the left (right) figure corresponds to M =P (As). We will come back to this point in the next subsection. The volume dependence of the DOS at the Fermi level, DOS(E F ), varies in accord with W SOMO , not W SH , as seen in figure 4(c). It is natural that DOS is inversely related to W SOMO because this directly measures the width of the band that intersects the Fermi level, but the point here is that W SOMO is not directly correlated with W SH and thus the unit cell volume. Figure 5 shows the effective tight-binding model adopted to fit the ab-initio band. In addition to the nearest-neighbor transfer energies (left part in figure 5, same notation is used as in previous studies), we need to introduce the next-nearest-neighbor transfer energies (right panel in figure 5) to reproduce the ab-initio band more accurately. Note that the stacking direction of the BDA-TTP molecules is taken in the c + a-direction of the lattice structure in our study, as shown in figure 1(b). We show in figure 6(a) the band structure of the two band tight-binding model for β-(BDA-TTP) 2 SbF 6 superposed to the ab-initio band structure. It can be seen that the tight-binding model nicely reproduces the ab-initio band structure. figure 6(b) shows the Fermi surface obtained from the model, which is essentially the same as the ab-initio result in figure 2(c).

Effective tight-binding model
Similar tight-binding fit is performed for all the anions, and the thereby determined transfer energies are summarized in Table 1. The bottom two lines of Table 1 represent the magnitude of the dimerization, t p1 /t p2 , in the stacking direction, and the relative transfer perpendicular to that direction (i.e. the strength of the two dimensionality), (t q1 + t q2 )/(t p1 + t p2 ). We can now quantitatively discuss the effect of the dimerization on the unit-cell volume dependence of W SOMO and DOS(E F ). Table 1 shows that t p1 /t p2 obtained by using the structural data at room temperature is around 1.22 except for M =As, in which t p1 /t p2 is 1.09. This means that the dimerization gap between the SOMO and HOMO bands is small in As, hence resulting in a larger W SOMO despite  the smaller W SH . This is consistent with the qualitative discussion given in subsection 3.1. From the viewpoint of the transfer energies, we consider that the difference of the anisotropy of the Fermi surface is due to the magnitude of the two-dimensionality, (t q1 + t q2 )/(t p1 + t p2 ). Our result shows that the average value of (t q1 + t q2 )/(t p1 + t p2 ) over the salts listed in Table 1 is around 0.46, which is smaller than that obtained by the extended Hückel calculation [3,4,15,17]. Therefore, the present Fermi surface exhibits an anisotropy in which the a + c-direction is more conductive than in the extended Hückel result.

Random phase approximation and the superconducting gap equation
We now move on to the RPA analysis for the effective models. To focus on the anion dependence, here we concentrate on the models obtained using the lattice parameters at room temperature. We show in figure7(a) the spin susceptibility in the Hubbard model of β-(BDA-TTP) 2 SbF 6 , where U = 0.38[eV] and the temperature T = 0.008[eV] here. The wave vector at which the spin susceptibility is maximized is k = (0.45π, 0.53π). The arrows in figure7(b) corresponds to the peak position of the spin susceptibility in (a), which can be considered as the nesting vector of the Fermi surface. In fact, there is a certain overlap between the original Fermi surface and the ones translated by these arrows, as shown in figure 7(d). The superconductivity is analysed by solving the in previous theoretical studies that adopt extended Hückel models (see the appendix) [15,16,17], and hence is not consistent with the STM study [10]. On the other hand, the present result is indeed more consistent with the recent measurement of the anisotropy of H c2 in β-(BDA-TTP) 2 SbF 6 [11]. Further study is necessary to resolve the consistency between existing experiments and theories. T c is obtained for the four anions for various on-site interaction U. figure 8 shows T c as functions of the unit-cell volume. There are several features found in the calculation results: (i) T c is the highest for Sb, (ii) T c somewhat decreases as we go from P to As, and (iii) T c goes down as we go from Sb to Ta, but it is still finite, close to that of P. These features are in accord with the anion dependence of DOS(E F ) or W SOMO shown in Fig 4(c) or (b), but not W SH or the unit-cell volume. Point (i) is indeed consistent with the experimental finding. Point (ii) can also give an understanding for the experimental observation that T c does not increase substantially as we go from P to As. This may at first sight be puzzling from the viewpoint of the increase in the unit-cell volume, from which one would expect an increase in the density of states. The key for the present understanding is the difference in the strength of the dimerization between P and As as mentioned previously. As for point (iii), T c of Ta being somewhat close to P, is not consistent with the experiments, where no superconductivity is observed for Ta. The discrepancy may be due to effects that are not taken into account in the present model (such as phonons), or in the RPA (more sophisticated electron correlation effects). On the other hand, the decrease of the T c with the increase of the unit-cell volume found in the present calculation might partially be related to the loss of superconductivity.

Conclusion
In the present paper, we have obtained the ab-initio band structure and the effective tight-binding model of β-(BDA-TTP) 2 MF 6 (M =P, As, Sb, Ta), and studied the superconducting gap form as well as the anion dependence of T c assuming spin fluctuation mediated pairing.
The results of ab-initio band calculation for β-(BDA-TTP) 2 SbF 6 show that the direction of the major axis of the cross section of the Fermi surface lies in the Γ-B direction, in contrast to the extended Hückel result. The density of states at the Fermi level is correlated with the band width of SOMO rather than SOMO+HOMO, and hence is not necessarily correlated with the unit-cell volume.
To construct an effective tight-binding model that accurately reproduces the abinitio band structure, we need both nearest and next nearest neighbor transfer energies. From the obtained transfer energies, the dimerization of M =As is smaller than that of M =P and Sb, and hence the SOMO band width of M =P is smaller than that of M =As despite the larger SOMO+HOMO band width.
The RPA calculation shows that the wave vector at which the spin susceptibility is maximized lies in the a + c direction. The gap function of the superconductivity mediated by the spin fluctuations changes sign along the Fermi surface like in a d-wave gap, and the hot spots of the gap (gap maxima) is in the k a ±k c direction, again different from the extended Hückel model. T c of M =Sb is the highest among all the anions in qualitative agreement with experiments.
There are some remaining issues that should be resolved in future studies. As mentioned above, the relation between the present Fermi surface, those obtained by extended Hückel calculation [3,4,15,17], and those obtained in two experiments [5,7] remains as an open question. As for the band calculation, we consider that the difference between the present ab-initio result and the extended Hückel result is intimately related to the empirical parameter in the extended Hückel calculation that controls including/excluding the 3d-orbital of the sulfur atoms of the TTP skeleton donors [22,23,24]. Actually, the extended Hückel result in which the 3d-orbital of the sulfur atoms is excluded [25] is similar to the present ab-initio result.
Although all experimental and theoretical studies suggest anisotropic superconducting gap, the form of the gap function is different between the present and previous theoretical studies based on extend Hückel band [15,17]. Given the consistency between the recent experiment [7] and the present calculation regarding the direction of the anisotropy of the Fermi surface, we believe that the gap maxima should lie in the k a ±k c direction as far as the spin-fluctuation-mediated pairing mechanism is concerned. If the nodes of the superconducting gap do not lie in this direction in the actual materials, then the present study would conversely suggest that a different pairing mechanism may be at work. In this sense, it is highly important to resolve the experimental controversy regarding the direction of the nodes in the superconducting gap [10,11].