Electronic Structure Calculation and Superconductivity in $\lambda$-(BETS)$_{2}$GaCl$_{4}$

Quasi-two-dimensional molecular conductor $\lambda$-(BETS)$_2$GaCl$_4$ shows superconductivity (SC) below 5.5K, neighboring the dimer-type Mott insulating phase. To elucidate the origin of SC and its gap function, we carry out first-principles band calculation and derive a four-band tight-binding model from the maximally localized Wannier orbitals. Considering the spin-fluctuation-mediated mechanism by adding the Hubbard $U$-term to the model, we analyze the SC gap function by applying the random phase approximation. We show that the SC gap changes its sign four times along the Fermi surface (FS) in the unfolded Brillouin zone, suggestive of a $d$-wave-like SC gap, which only has two-fold symmetry because of the low symmetry of the crystal structure. Decomposing the SC gap into the pairing functions along the crystal axes, we compare the result to similar analysis of the well-studied $\kappa$-type molecular conductors and to the experiments.

The quasi-two-dimensional (Q2D) molecular conductor λ-(BETS) 2 GaCl 4 , where BETS is bis(ethylenedithio)tetraselenafulvalene, exhibits superconductivity (SC) below 5.5 K. [1][2][3] It has attracted interest as a candidate for realizing the FFLO state under magnetic field owing to its highly two-dimensional electronic structure. [4][5][6] Another interest is that its isostructural compound λ-(BETS) 2 FeCl 4 shows a field-induced SC phase under strong magnetic field, 7,8) considered to be connected to that of the Ga salt (at zero field) as indicated by the measurements of alloyed samples. 9) Despite extensive experimental works, theoretical investigation of SC in this compound from a microscopic viewpoint has been lacking, which is the purpose of this study.
In the λ-type structure, the BETS molecules stack along the a direction, forming a triclinic unit cell with space group P1. 1,2) There are dimers of BETS molecules with large intradimer transfer integrals (termed t A ), which show further dimerization, i.e., a tetramer of BETS forms the unit cell. The GaCl −1 4 closed-shell anion sheets lead to the highest-occupied molecular orbital (HOMO) of BETS forming a Q2D quarter-filled system in terms of holes. 3) From its dimerized structure, whose limit of large dimerization will be a half-filled system, the electronic structure has an analogy with the well-studied Mott transition system κ-(ET) 2 X [ET = bis(ethylenedithio)tetrathiafulvalene]. In fact, by chemical substitution in the anions GaX z Y 4−z (X, Y =F, Cl, Br) 10,11) or by choosing different donor molecules, 12) the SC phase is suggested to locate next to the Mott insulating phase as in κ-(ET) 2 X. 13) Although the nature of the insulating state just near the SC phase remains to be clarified, i.e., whether it is * aizawa@kanagawa-u.ac.jp non-magnetic 10,14) or antiferromagnetic, 12,15) a recent NMR measurement in λ-(BETS) 2 GaCl 4 reports the development of spin fluctuations above the SC transition temperature. 16) As for the SC gap function, early measurements show a two-fold symmetry within the conductive plane, by means of the anisotropy of the upper critical field H c2 17) and of the flux-flow resistivity. 18) Ref. 18 observes a dip structure in the angle dependence of the resistivity under magnetic field, when the magnetic field is applied parallel to the c axis. More recently, a heat capacity measurement indicated the line-nodal gap of dwave pairing, 19) whereas a µSR measurement reports a possible mixture of the extended s-and d-wave gaps. 20) The electronic structure of λ-(BETS) 2 GaCl 4 has been discussed within the tight-binding model based on the HOMO of the BETS molecule, where the transfer integrals are calculated using the extended Hückel method. 11,21,22) The band structure near the Fermi energy shows four bands since the unit cell contains four BETS molecules as mentioned above. The calculated Fermi surface (FS) is similar to that of the κ-(ET) 2 X, 23) which consists of a pair of open and closed FS, despite the difference in their molecular packings. One issue is that, since the extended Hückel method contains semiempirical parameters, there are estimates with appreciable discrepancies.
In this study, we present the band structure obtained from first-principles calculations, and estimate the transfer integrals of the four-band model from the maximally localized Wannier orbitals (MLWO). Then, considering the pairing mechanism mediated by the spin fluctuations, we apply the random phase approximation (RPA) to the four-band Hubbard model of λ-(BETS) 2 GaCl 4 . The results show a d-wave-like SC gap and we will discuss its origin related to the spin susceptibility.
The first-principles band calculations were performed within density functional theory (DFT) with generalized gradient approximation 24) using WIEN2k, 25) and a tight-binding model was derived by applying MLWO 26,27) scheme using wannier90 package. 28) Figure 1(a) shows the two-dimensional band dispersions near the Fermi level for the experimental structure data. 11) Dispersion along the interlayer direction is small, of the order of 0.1 meV. There are four bands originated from HOMO of BETS, corresponding to the extended Hückel bands. One point we note is that, since nearly flat band dispersions are present near Z point, the density of states (DOS) exhibits a van-Hove singularity (vHS) slightly below the Fermi level, as shown in Fig. 1(b). In Fig. 1(c), we show the FS obtained from the DFT calculation, consisting of open and closed portions, which we call FS0 and FS1 in the following. The former comes from the top band and the latter comes from the second-totop band. The shape of the FS is similar to the extended Hückel results. 11,21,22) We regard these four bands as the target bands and derive a tight-binding model by constructing a MLWO on each molecule. As shown in Fig. 1(a), the band structure of the four-band model, which includes the distant transfer integrals, accurately reproduces the DFT band dispersion. Table I. Transfer integrals and site-energy difference in meV for λ-(BETS) 2 GaCl 4 , where the site-energy difference between the BETS-1(4) and BETS-2(3) is defined as ∆E ≡ E 2(3) − E 1 (4) . The superscript, eH, stands for the extended Hückel results 11,21,22) and the subscript, Fe, represents the results of λ-(BETS) 2 FeCl 4 , having the same crystal structure. We summarize the obtained transfer integrals and the energy difference between the nonequivalent BETS, ∆E, in Table I, together with the extended Hückel results in the literature. The notation of inter-molecular bonds is as shown in Fig. 1(d). The other transfer integrals not listed here have absolute values less than 13 meV. In the RPA analysis below, we use all the obtained transfer integrals in the two-dimensional plane. 29) As a common feature among our results and the previous extended Hückel calculations, t A is the largest, which is the intradimer transfer integral. This gives the splitting between the upper two and lower two bands, approximately corresponding to the antibonding and bonding combinations of the HOMO. The transfer integrals along the stacking direction t B and t C have close values in contrast with previous data, which indicates that the degree of tetramerization is smaller than previously discussed. 14) The effective transfer integrals between the antibonding combination of HOMO of BETS dimers along the a direction can be approximated from the large dimerization limit ast B ≡ t B /2 andt C ≡ t C /2; that in the transverse direction is t ⊥ ≡ (t p + t q + t r ) /2. Our results show a relation |t B | ≃ |t C | ≃ t ⊥ . Then, the BETS dimers possess a square-lattice-like network along the a and c directions, with weaker diagonal transfer integrals We can interpret the large DOS to be originated from this relation since the ideal square lattice has a singularity of the DOS at half-filling. Another recent DFT calculation based on the pseudopotential shows the same result. 20) Next, by introducing the on-site Coulomb interaction U to the four-band model, we study the spin susceptibility χ sp and the SC gap function within the framework of the spin-fluctuation-mediated pairing mechanism. The Hamiltonian is described as where i and j are unit-cell indices, α and β specify the sites 1-4 in a unit cell [see Fig. 1 (d)], c † iασ (c iασ ) is the creation (annihilation) operator for spin σ at site α in unit cell i. t iα:jβ is the transfer integral between site (i, α) and site (j, β), estimated as above, and iα : jβ represents the site pairs. n iασ is the number operator for electrons with spin σ on site α in unit cell i and n iα = n iα↑ + n iα↓ .
To deal with the effect of the Coulomb interaction U , we apply the multisite RPA, e.g., described in Ref. 30; here we focus on SC state in a situation where other instabilities are weaker. The Green's function, as well as the susceptibilities, pairing interaction, and SC gap function are all 4×4 matrices. The gap functionφ (k , iε n ) and its eigenvalue λ are obtained by solving the linearized Eliashberg equation. The critical temperature T c corresponds to the temperature where λ reaches unity. Because we consider only the on-site interaction U , the spin susceptibility is much larger than the charge susceptibility. Therefore, we will show the spin susceptibility χ sp obtained from the largest eigenvalue for the lowest Matsubara frequency. The SC gap function is presented in the band representation at the lowest Matsubara frequency. In the present calculation, we take 96×96 k-point meshes and 16384 Matsubara frequencies. The on-site interaction is chosen as U = 0.4 eV.
As shown in Fig. 2(a), at temperature T = 0.006 eV (≃ 70 K), the spin susceptibility χ sp has the maximum value around Q 0 = (Q 0a , Q 0c ) ≃ (−3π/8, 3π/8) and a broad substructure around Q 1 = (Q 1a , Q 1c ) ≃ (−π/6, 5π/6). To discuss the nesting properties in the following, χ sp in the extended zone along the Γ-Z direction is shown and we defineQ 0 = −Q 0 + (0, 2π). In Fig. 2(b), we show the FS on the left and the SC gap function for the top (second-to-top) band, namely, for FS0 (FS1) on the center (right) for λ ≃ 0.42 31) . The FS can approximately be regarded as an ellipse in the unfolded Brillouin zone, whereas FS0 and FS1 are slightly disconnected around the k points where they approach to each other. In the following, we will call this point as the crossing point and the elliptic FS in the unfolded zone as the 'extended FS'. As we can see in the figure, wave-number vectors Q 0 and Q 1 correspond to the FS nesting. As for the SC gap function, first we note that for FS0 (FS1) it has a positive (negative) sign along almost the whole Brillouin zone. This gives rise to the large s-wave components in the analysis below.
To clarify the relation between the electronic structure and the SC gap, in Fig. 2(c) we plot the SC gap functions for k vectors within ±0.01 eV from the Fermi level, which corresponds to be about 1% of the band width (of four bands), from the Fermi level in the extended zone. Then, we can see that the SC gap changes its sign four times along the extended FS reminiscent of a d-wavelike gap, which possesses two kinds of nodes, from which gentle/steep increase of the SC gap is seen. We call them as gentle/steep node structures. In the figure, the two nesting vectorsQ 0 and Q 1 connect the portions of FS where the SC gap has a large amplitude and shows sign changes. The positions of the FS where the SC gap amplitude is large almost coincide with the positions giving rise to the vHSs, resulting in a high stability of the gap structure. In fact, a similar analysis based on the effective two-band model (the upper two bands) also gives rise to similar angular dependence. 20)  Next, we attempt to decompose the SC gap into different components, as has been done for κ-(ET) 2 X. [32][33][34][35] In the λ-type structure, the point group symmetry is low as C i , therefore naturally different components mix. 36) Therefore, we decompose the d-wave-like gap into pairing components along the crystal axes. Here we take crystal c and a directions as x and y axes, respectively, to make the correspondence between other systems clearer. In this choice of axes, the SC gap structure in Fig. 2 (c) apparently looks close to a d x 2 −y 2 -wave gap since the nodes are nearly along the diagonal directions. We introduce the fitting function given as c+a cos(k x + k y ) + C f c−a cos(k x − k y ) + · · · (up to 20th nearest neighbors), (2) where f is for the choice of the two bands represented by FS0 or FS1, the subscript represents the pairing direction, and the fitting variables C f are the weights of the basis function on the FS "f ". Longer range pairing states as represented in Eq. (2) are also considered in the actual calculation, but have small contributions.
We summarize the ratio of the fitting variables for the basis function in Table II. In the case of FS0 [center of Fig. 2 (b)], although the ratio in the c − a direction is the largest, the pairing ratio along the a direction is subdominant and comparable with that of the intra-unit cell. It is suggestive that the SC gap of the FS0 is affected by the pairing along the a direction, in which the BETS molecules stack. By contrast, for FS1, the three components, namely the intra-unit cell as well as c and a directions, are comparable. As expected, the SC gap of FS1 exhibits a two-dimensional pairing. Table II. Ratio of the fitting variables of the basis function on the FS "f ", from the d-wave-like gap function. We take the intraunit-cell component C f 0 as unity; to stress the different sign between the two bands, we put different signs.

Fitting variable
FS0 FS1 To compare with the previous studies of κ-(ET) 2 X, [32][33][34][35] we rewrite the ratio of the well-known SC gap, as d x 2 −y 2 -, d xy -, extended s 1(2) -wave, which is the pairing with the same sign between the first (second) nearest neighbors. Note that we decompose the d-wavelike gap into the well-known SC gap and confirm that the same components are obtained. We list the components of the SC gap in Table III. Several SC-gap components of the FS0 are comparable. By contrast, for the FS1, the components of the isotropic s-and extended s 1 -wave possess large negative value. We should note that, even though the component of the isotropic s-wave, which is same as the intra-unit-cell pairing, is large, this does not mean that the pairing, in real space picture, occurs on the same BETS molecule, since the SC components are obtained in the "folded" Brillouin zone. Namely, an anisotropic pairing, e.g., the nearest neighbor pairing between BETS-2 and BETS-1 or BETS-3 within the same unit cell, is converted to an isotropic s-wave component in the folded Brillouin zone because the pairing occurs within the unit cell. The results here that multiple components have comparable values are noticeably different from the case of κ-(ET) 2 X. [32][33][34][35] In that case, the effective half-filled dimer Hubbard model shows the instability toward d x 2 −y 2wave SC in the extended zone 32,[37][38][39][40][41][42] while for the 3/4filled model realistic parameters provide d xy -type -wave gap [33][34][35][41][42][43] but with considerable extended s-wave component. [32][33][34][35] We can attribute such difference to the different crystal structure geometries: κ-type has a D 2h point group symmetry, so that pure d x 2 −y 2 -wave can be stabilized but not pure d xy -wave in the extended zone. κ-(ET) 2 X has parameters close to the triangular lattice giving rise to geometrical frustration effect, while our analysis here provides a square-lattice like network, as in the high T c cuprates producing the stability of d x 2 −y 2wave, but with large mixing with other components of the well-known SC gap due to the low symmetry of the crystal structure.
Finally let us discuss the experimental works from the viewpoint of our results giving the d-wave-like gap. The results in the transport measurements indicating the twofold symmetry of the angular dependence of the SC gap in this compound are compatible with our results since the d-wave-like gap only possesses the two-fold symmetry. 17,18) The existence of the nodal SC gap is suggested by a recent measurement of the heat capacity, 19) which is in accordance with our results showing nodes along the diagonal directions. As for the nodal position, the fluxflow resistivity measurement suggests a dip structure of the resistivity when the magnetic field is applied parallel to the c axis. 18) This is consistent with the d-wave-like gap we obtained, namely, the large SC gap around vHS and the steep node structure are present. A recent µSR measurement suggests that the SC of this compound is a mixture of the extended s-wave and d-wave SC. 20) A direct comparison between our results might be difficult since the method of decomposing the SC gap is different from ours here, magebut nevertheless, the mixture of different components is indeed consistent.
In conclusion, we have obtained the DFT band structure and the four-band model of the Q2D molecular conductor λ-(BETS) 2 GaCl 4 . Within the spin-fluctuationmediated pairing mechanism, we study the SC gap function and its properties by applying the RPA. The network of the BETS dimers shows a square-lattice-like structure, giving rise to large DOS near the Fermi level. We propose that the FS nesting within this characteristic electronic structure results in the d-wave-like SC gap, which changes its sign four times along the extended FS and possesses the two-fold symmetry.
To elucidate the pairing components of the d-wave-like gap, we have decomposed this gap function into the pairing components along the crystal axes, and estimate the pairing ratio for each FS. We have shown that the SC gap of FS0 is affected by the pairing in the stacking direction of the BETS and that the gap of FS1 exhibits a twodimensional pairing. To compare the previous studies of κ-(ET) 2 X, we transform the component of the pairing in the crystal axes to that of the well-known SC gap functions, and show that the several SC gap components can be comparable in both FSs.
The effect of strong electronic correlation beyond RPA, which is expected to play a role since the system is con-sidered to be located near the Mott transition, is an interesting issue left for future studies.