$d$-wave Superconductivity in the Hubbard model on the isotropic triangular lattice and a possibility of the chiral $d+id$ pairing as a quasi-stable state

Strong electron correlations lead to interesting phenomena such as superconductivity with various pairing symmetries and purely paramagnetic insulator (spin liquid) in low dimensional materials. The organic charge-transfer salts $\kappa$-(BEDT-TTF)₂Cu₂(CN)₃ [1, 2, 3, 4, 5, 6, 7] is a good example of such materials, whose spin liquid state transits to a superconductor with $d_{xy}$ pairing symmetry at $T_{K}=3.9$ K upon applying pressure [5, 6].

The superconductivity(SC) of this material was studied theoretically using the Hubbard model on the isotropic triangular lattice, which is a simple effective Hamiltonian of this material[8] described by the hopping $t$ and the on-site Coulomb repulsion $U$. The variational Monte Carlo[9], exact solution[10], renormalization group[11, 12], and density renormalization group[13] were applied to this model, and these studies excluded the $d$-wave SC below the Mott transition point at half-filling. These results imply that the SC of this material can not be described by the Hubbard model on the triangular lattice.

Contrary to these analyses, the studies[14, 15, 16] using the cluster dynamical mean filed theories(CDMFT) predicted the $d$-wave SC. In the study of the cellular dynamical mean-field theory[14], 2$\times$2 cluster is used as the reference cluster and 8 bath sites are attached to that, and it is found that the $d_{x^{2}-y^{2}}$ SC is realized below the Mott transition point. The analysis of the variational cluster approximation(VCA)[15] adopted 2$\times$2 and 2$\times$4 clusters as the reference clusters and reported that the $d_{x^{2}-y^{2}}$ SC is the ground state near the Mott insulator phase and the gap symmetry changes to $d_{xy}$ for lower values of $U$. However the time-reversal symmetry breaking chiral $d+id$ SC, which is an important candidate of the pairing symmetry on the triangular lattice, was not considered in these analyses. The chiral $d+id$ state was investigated by VCA using 2$\times$2 cluster and 6-site triangular clusters[16] as the reference clusters, and it is found that the 2$\times$2 site analysis yields a strong preference for the $d_{x^{2}-y^{2}}$ SC, while the ground state is the chiral $d+id$ SC below the Mott transition point on the 6-site triangular cluster. So the pairing symmetries predicted in the CDMFT do not agree with the experiments[5, 6], which suggests together with the results of the anayses[9, 10, 11, 12, 13], that some physics factors not included in the simple Hubbard model may be necessary to understand the SC of $\kappa$-(BEDT-TTF)₂Cu₂(CN)₃ and related materials. Even within the analyses of the CDMFT, the pairing symmetry of the SC phase is still controversial.

However, in the CDMFT like VCA, it is guaranteed that the results converge to the exact results as the size of the reference cluster increases because the electron correlations within the reference cluster are exactly taken into account. Therefore we improve the previous studies of the CDMFT[14, 15, 16] by increasing the reference cluster size at least until the semi-quantitative convergence is observed between two different sizes of the reference clusters. After this minimum check, we compare the obtained results with the experiments[5, 6] to see if the SC of this material is described by this Hamiltonian.

In this paper we study the $d$-wave SC in the Hubbard model on the isotropic triangular lattice by VCA using the 12-site and 14-site clusters in Fig. 1 at zero temperature and half-filling. We found that the results are semi-quantitatively the same for 12-site and 14-site clusters and the ground state is the $d_{xy}$ SC below the Mott insulator phase. The energy difference between the normal paramagnetic state(PM) and the $d_{xy}$ SC is $0.02t\sim 0.04t$ for $U/t\sim 5$. The energy of the chiral $d+id$ state is about $0.01t$ higher than that of the $d_{xy}$-wave ground state for $U/t\sim 5$. The $d$-wave SC is not realized above the Mott transition point $7\lesssim U/t$. We have confirmed that our overall conclusion is not changed by slight differences of the probes of the symmetry breaking patterns by using the superconducting form factors (e) and (f) in addition to (b)$\sim$(d) in Fig. 1 into our analysis.

Comparing our results with the experiments of $\kappa$-(BEDT-TTF)₂Cu₂(CN)₃, our prediction of the gap symmetry $d_{xy}$ coincides with the analysis of the thermal conductivity[5] and the STM observations[6]. Adopting the estimate $t\sim 0.06$ eV[17, 18, 19], the energy difference $0.02t\sim 0.04t$ between the PM and $d_{xy}$ SC is semi-quantitatively consistent with the transition temperature $T_{K}=3.9$ K, as will be discussed in detail later. Thus we consider that the SC of this material is well-described by this model. Because of the semi-quantitative convergence of our results and agreement with the experiments, we expect that our overall conclusions in VCA are robust with respect to a further increase in the cluster size.

Our prediction of the ground state pairing symmetry disagrees with the preceding analyses of the CDMFT[14, 15, 16] with smaller reference clusters and we shall discuss in detail about the origin of the discrepancies later, considering a relation of the shapes of reference clusters and symmetry of the system using the concept of the effective potential.

The Hamiltonian of the Hubbard model on the isotropic triangular lattice reads

where $t_{ij}=t$ for the solid lines in Fig. 1, $U$ is the on-site Coulomb repulsion, and $\mu$ is the chemical potential. The annihilation (creation) operator for an electron at site $i$ with spin $\sigma$ is denoted as $c_{j\sigma}$ ($c_{i\sigma}^{\dagger}$) and $n_{i\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}$. The energy unit is set as $t=1$ hereafter.

We use VCA[20, 21, 22] in our analysis. In this approach we start with the thermodynamic grand-potential $\Omega_{\mathbf{t}}$ written in the form of a functional of the self-energy $\Sigma$ as

In Eq. (2), $G_{0}$ is the non-interacting Green’s function of $H$, $F[\Sigma]$ is the Legendre transform of the Luttinger-Ward functional[23], and the index $\mathbf{t}$ denotes the explicit dependence of $\Omega_{\mathbf{t}}$ on all the one-body operators in $H$. The variational principle $\delta\Omega_{\mathbf{t}}[\Sigma]/\delta\Sigma=0$ leads to the Dyson’s equation. Eq. (2) gives the exact grand potential of $H$ for the exact self-energy of $H$, which satisfies Dyson’s equation.

All Hamiltonians with the same interaction part share the same functional form of $F[\Sigma]$, and using that property we can evaluate $F[\Sigma]$ for $\Sigma$ of a simpler Hamiltonian $H^{\prime}$ by exactly solving it. In VCA, we divide the original infinite cluster into the identical clusters, referred to as the reference cluster, that tile the infinite lattice, and construct $H^{\prime}$ from $H$ by removing the hopping interactions between these clusters. Then writing the Eq. (2) for $H^{\prime}$, and substituting it from $\Omega_{\mathbf{t}}[\Sigma]$ of $H$, we obtain

where $G^{\prime}_{0}$ is the non-interacting Green’s function of $H^{\prime}$ and $\mathbf{t}^{\prime}$ denotes all the one-body operators in $H^{\prime}$. In Eq. (3) we evaluate $\Omega^{\prime}_{\mathbf{t^{\prime}}}[\Sigma]$ for $\Sigma$ of $H^{\prime}$ by exactly solving it, then $\Omega_{\mathbf{t}}[\Sigma]$ becomes a function of $\mathbf{t}^{\prime}$ expressed as

where $\Omega^{\prime}_{\mathbf{t^{\prime}}}$ is the exact grand potential of $H^{\prime}$ and the functional trace has become an integral over the diagonal variables (frequency and super-lattice wave vectors) of the logarithm of the determinant over intra-cluster indices. The frequency integral is carried along the imaginary axis and $\delta\rightarrow+0$.

The variational principle $\delta\Omega_{\mathbf{t}}[\Sigma]/\delta\Sigma=0$ is reduced to the stationary condition $\delta\Omega_{\mathbf{t}}(\mathbf{t}^{\prime})/\delta\mathbf{t}^{\prime}=0$, and its solution and the exact self-energy of $H^{\prime}$ at the stationary point, denoted as $\Sigma^{*}$, are the approximate grand-potential and self-energy of $H$ in VCA. Physical quantities, such as expectation values of one-body operators, are evaluated using the Green’s function $G_{0}{}^{-1}-\Sigma^{*}$. In VCA, the restriction of the space of the self-energies $\Sigma$ into that of $H^{\prime}$ is the only approximation involved and short-range correlations within the reference cluster are exactly taken into account by exactly solving $H^{\prime}$. A possible symmetry breaking is investigated by including in $H^{\prime}$ the corresponding Weiss field that will be determined by minimizing the grand-potential $\Omega_{\mathbf{t}}$.

In our analysis, the 12-site and 14-site clusters in Fig. 1 are used as the reference clusters to set up the cluster Hamiltonian $H^{\prime}$. We refer to these clusters as 12D and 14D hereafter. Within these reference clusters, every site is connected to at least three other sites.

To study the superconductivity, we include the Weiss field Hamiltonian

into $H^{\prime}$, where $z$ is treated as a variational parameter, and adopt the Nambu formalism $c_{j\uparrow}=\tilde{c}_{j\uparrow}$ and $c_{i\downarrow}=\tilde{c}_{i\downarrow}^{\dagger}$.

We classify the real-space superconducting form factors $\Delta_{ij}$ between the nearest-neighbor sites into the irreducible representations of the invariant group of regular hexagon $C_{6v}$, and take as (b) $d_{1}$ and (c) $d_{2}$ in Fig. 1. $C_{6v}$ consists of six rotations (one is identity) and six mirror reflections, and $d_{1}$ is anti-symmetric under two of the six mirror reflections while $d_{2}$ is symmetric under these reflections. As for the time-reversal symmetry breaking chiral superconducting state, we take the combination (d) $i(\sqrt{3}d_{1}+id_{2})/2$ and denote as $d_{1}+id_{2}$ hereafter.

In the STM experiments[6] the two line nodes of the gap are ${\pi}/{4}$ from the axes of the two reciprocal lattice vectors, which indicates $d_{2}$, whose two line nodes lie between the basic lattice vectors in the real space. Also, the observed gap[6] is symmetric under two mirror reflections. In the case of $d_{1}$, one of the two line nodes lies along the basic lattice vectors, while the other is perpendicular to that, which does not agree with the result of STM.

To confirm that subtle differences of the choice of the Weiss field do not change our results, we also consider (e) $d^{\prime}_{2}$ in Fig. 1, which involves next-nearest neighbor interactions, and (f) $d_{1}+id^{\prime}_{2}$. $d^{\prime}_{2}$ has the same symmetric property as $d_{2}$ and also probes the SC in the $d_{2}$ direction. The two chiral states $d_{1}+id_{2}$ and $d_{1}+id^{\prime}_{2}$ are topologically equivalent since both the real and imaginary parts change the sign twice in a similar manner by $2\pi$ rotation. Because the terminology of the gap $d_{1}\equiv d_{x^{2}-y^{2}}$ and $d^{\prime}_{2}\equiv d_{xy}$ on the square lattice are widely used in the previous studies[5, 6, 14, 15, 16], we also denote $d_{1}$ as $d_{x^{2}-y^{2}}$, and both $d_{2}$ and $d^{\prime}_{2}$ as $d_{xy}$.

With this set up, $\mathbf{t}^{\prime}$ is reduced to the Weiss parameter $z$ in Eq. (5) and the cluster chemical potential $\mu^{\prime}$, where $\mu^{\prime}$ should be included for the thermodynamic consistency[24], and the stationary condition $\delta\Omega_{\mathbf{t}}(\mathbf{t}^{\prime})/\delta\mathbf{t}^{\prime}=0$, is solved by searching the stationary point of $\Omega(\mu^{\prime},z)$, which we denote as the grand-potential per site. During the search, the chemical potential of the system $\mu$ is also adjusted so that the electron density $n$ is equal to 1 within the accuracy of $10^{-5}$. The energy per site is given by $E=\Omega+\mu n$ where $\Omega$ is the value of $\Omega(\mu^{\prime},z)$ at the stationary point. In general, a stationary solution with $z\neq 0$ corresponding to the superconducting state and that with $z=0$ corresponding to the normal paramagnetic state are obtained, and these energies are compared to determine the ground state.

Before the analysis of the SC, we compute the Mott transition point for the PM. To examine the gap we calculated the density of state per site

imposing $z=0$, where the $k$ integration is over the corresponding Brillouin zone, and $\eta\rightarrow 0$ limit is evaluated using the standard extrapolation technique, and obtained that the Mott transition point $U_{M}$ is $U_{M}=6.3$ on 12D[16, 25, 26], and $U_{M}=5.4$ on 14D.

Next we consider the SC for $U_{M}\lesssim U$. At $U=9$ and $U=12$, we found that $\Omega(\mu^{\prime},z)$ is monotonically increasing as a function of $z$ near the half-filling and no SC solutions are obtained. At $U=7$, we found the stationary solution for all $d_{1}$, $d_{2}$, $d_{1}+id_{2}$, $d_{1}+id^{\prime}_{2}$, and $d^{\prime}_{2}$ on both 12D and 14D, and their energies are degenerate with the PM within about $10^{-4}t\sim 10^{-5}t$. Therefore we consider that the SC is not realized for $U_{M}\lesssim U$.

Next we consider the SC for $U\lesssim U_{M}$. Figs. 2–4 show $\Omega(\mu^{\prime},z)$ as functions of $z$ and $\mu^{\prime}$ computed by VCA on 14D at $U=4$ for $d_{2}$, ${d^{\prime}}_{2}$, $d_{1}+i{d}_{2}$ and $d_{1}+i{d^{\prime}}_{2}$. In these figures, the marks show the stationary point satisfying the half-filling condition $n=1$. In figures (a), we fix $\mu^{\prime}$ to be the value of the stationary solution, and in (b) we fix $z$ to be the value of the stationary solution. The values of $\mu$ are given in the caption. Similarly, we found the stationary solutions satisfying $n=1$ for $d_{2}$, $d_{1}+id_{2}$, and $d_{1}+id^{\prime}_{2}$ at $U=1,2,3,4,5$,and $6$ on 12D and 14D. As for $d^{\prime}_{2}$, we found the stationary solutions for $1\leq U\leq 6$ on 14D and $4\leq U\leq 6$ on 12D, but we were not able to obtain the stationary solutions for $U\leq 3$ on 12D because $\Omega(\mu^{\prime},z)$ becomes discontinuous. In general, long range correlations become more important for smaller $U$, and 12D cluster is slightly small to simulate $d^{\prime}_{2}$ for $U\leq 3$. On 12D we continued the stationary point search up to larger values $z\simeq 1.6$. Sometimes these searches are terminated by the discontinuity of the grand potential due to the change of the electron numbers of the cluster ground state. By these searches we confirmed that the obtained solutions are unique.

As for $d_{1}$, we found that $\Omega(\mu^{\prime},z)$ remains monotonically decreasing as we increase $z$, until $\Omega(\mu^{\prime},z)$ becomes discontinuous both on 12D and 14D for $1\leq U\leq 6$. Therefore the stationary solution was not obtained for $d_{1}$. A monotonically decreasing behavior of the grand potential is also observed e.g., for $s$-wave and $d_{xy}$ pairings for the superconductivity in the Hubbard model on the square lattice in VCA using $2\times 2$ cluster[27]. Since we will be able to exclude the quasi-stable $s$-wave SC in the Hubbard model based on the physics ground, the monotonically decreasing behavior does not necessarily imply that the reference cluster is too small to simulate the corresponding ordered state. In our case, we consider that $d_{1}$ is simulated well on 12D and 14D because the other four pairings are simulated well on these clusters, and the absence of the stationary solution means $d_{1}$ is not realized for $1\leq U\leq 6$. Rigorously speaking, whenever stationary solutions are not obtained, there always remains a possibility that the reference cluster is not suited to simulate the corresponding states.

Next we analyze the ground state energies of the SC. Fig. 5 shows the difference of the energy per site between the SC and PM, $\Delta E=E_{SC}-E_{PM}$ as a function of $U$ computed in VCA on (a) 12D and (b) 14D. The dotted lines indicate the Mott transition point $U_{M}$. At $U=5$ on 14D, $\Omega(\mu^{\prime},z)$ becomes almost flat in the $z$ direction near the stationary point for $d_{2}$ and the determination of the minimum might contain numerical errors, so we analyzed the energy differences also for $U=4.8$ and $U=5.2$. Taking into account the fact that $d_{2}$ and $d^{\prime}_{2}$ share the same symmetry property and are topologically equivalent, the general features are the same for 12D and 14D, though the cluster size dependence is not negligible for the magnitudes of the energy differences. For $1\leq U\leq 6$, $d_{xy}$ ($d_{2}$ and $d^{\prime}_{2}$) SC is the ground state and the chiral states $d_{1}+id_{2}$ and $d_{1}+id^{\prime}_{2}$ are energetically disfavored compared to $d_{xy}$. The energy difference between the PM and $d_{xy}$ SC is about $0.02t\sim 0.06t$ for $U\simeq U_{M}$, and decrease as $U$ decreases. The energy difference between the $d+id$ and $d_{xy}$ is about $0.01t\sim 0.03t$ for $U/t\simeq 5$

We compare our results with experiments of $\kappa$-(BEDT-TTF)₂Cu₂(CN)₃, which exhibits a transition from spin liquid to the $d_{xy}$ SC at $T_{K}=3.9$ K upon applying pressure[1, 2, 3, 4, 5, 6, 7]. Our prediction of the pairing symmetry is $d_{xy}$ and agrees with the experiments. For more quantitative comparisons, we assume that applying the pressure has the effect of increasing $t$ without affecting $U$, thus decreasing $U/t$. Then the Mott insulator region $5\lesssim U/t\lesssim 7\sim 8$ slightly above the SC is the candidate for the spin liquid, so we adopt the estimate $5\lesssim U/t\lesssim 8$[17, 18, 19]. Our analysis excludes the estimate $12\lesssim U/t\lesssim 15$[19] since the pressure in the experiments will not be able to push this region of $U/t$ down to the $d_{xy}$-wave SC phase $U/t\lesssim 7$.

The energy difference between the PM and $d_{xy}$ SC is about $0.02t\sim 0.06$ for $U\simeq U_{M}$. Even though the precise determination of the transition temperature requires the analysis of the entropy factor and finite temperature free energy, the transition temperature will have to be sizably lower than this energy difference to realize this $d_{xy}$ SC, since the entropy factor will be larger for the higher energy PM compared to the SC, and tend to increase the probability of the appearance of the PM. Assuming $t\simeq 0.06$ eV[17], this energy difference corresponds to $14\sim 42$ K, and the transition temperature $T_{K}=3.9$ K in the experiment is sizably lower than this value.

There exists fully-gapped $d+id$ slightly higher than the nodal gap $d_{xy}$, and might be observed as a quasi-stable state. It is an interesting challenge to observe the transition from $d+id$ to $d_{xy}$, or some effects of $d+id$ in $d_{xy}$ phase in experiments of this material.

Finally we discuss about the origin of the discrepancies between the preceding analyses[14, 15, 16] and ours. In general it is important in the CDMFT with relatively small reference clusters to choose appropriate shapes of reference clusters taking into account the symmetry of the system because the electron correlations within the reference clusters are exactly taken into account. In the case of isotropic triangular lattice, there is $C_{6v}$ symmetry, and the order parameters of $d_{x^{2}-y^{2}}$ and $d_{xy}$ transform according to the two dimensional irreducible representation $E_{2}$ of $C_{6v}$[28]. Because of this property, they mix in the singlet effective potential classified by the group symmetry, which leads to the possibility of the formation of $d+id$. Also, the two states $d_{x^{2}-y^{2}}$ and $d_{xy}$ are degenerate up to the fourth-order expansion of the effective potential in terms of the order parameters, and sixth-order terms are necessary to solve this degeneracy[29]. The electron correlations within the square shape 2$\times$2 and 2$\times$4 clusters used in the preceding studies[14, 15, 16] would not be able to create well the effective potential of this symmetry property and these analyses could not predict accurately the stable SC pairing state on the isotropic triangular lattice. The 6-site triangular cluster used in Ref. 16 and our 12D do not keep $C_{6v}$ symmetry, but they keep $C_{3v}$ symmetry, and the situations are the same to the case of $C_{6v}$ in the sense that the order parameters of $d_{x^{2}-y^{2}}$ and $d_{xy}$ transform according to the two dimensional irreducible representation $E$ of $C_{3v}$[28], and they mix in the singlet effective potential, leading to the possibility of the formation of $d+id$.

The origin of the discrepancies between the results of the 6-site triangular cluster in Ref. 16 and ours will be the difference of the cluster size. Our results suggest that, even though electron correlations within 4-site or 8-site clusters may be adequate to predict some physics quantities in VCA, clusters of more than ten sites are necessary to predict delicate physics property like pairing symmetry of the SC on the isotropic triangular lattice. As for the cluster size dependence, we have checked that our results of 12D are semi-quantitatively the same to those of 14D. Our 14D does not keep $C_{3v}$ symmetry, but the results of VCA converge to the exact results as the size of the reference cluster increases regardless of the detailed shapes of the cluster, and in fact they did for 12D and 14D. Since the physics discussed here seems to be delicate, we have also checked that our overall conclusions are not changed by the slight difference of the choice of the Weiss fields. Because of these checks and the semi-quantitative agreement between the experiments and our results, we expect that our results are semi-quantitatively robust with respect to a further increase in the reference cluster size.

We have studied the $d$-wave superconductivity in the Hubbard model on the isotropic triangular lattice by VCA at $n=1$. We have improved the preceding analyses by increasing the reference cluster size while paying attention to the cluster shape and symmetry of the system. We found that the results are the semi-quantitatively the same for the 12-site and 14-site clusters and the ground state is the $d_{xy}$ SC below the Mott transition point $U/t\lesssim 7$. The $d$-wave SC is not realized for $7\lesssim U$. Our prediction of the pairing symmetry $d_{xy}$ agrees with the electronic thermal conductivity[5] and STM experiments[6] of $\kappa$-(BEDT-TTF)₂Cu₂(CN)₃, and adopting the estimate $5\lesssim U/t\lesssim 8$ and $t=0.06$ eV[17], our result is semi-quantitatively consistent with the SC transition temperature $T_{K}=3.9$ K of this material. We also found that there exists the chiral $d+id$ SC slightly above the ground state $d_{xy}$ SC, and it is an interesting challenge to observe the transition from $d+id$ to $d_{xy}$, or some effects of $d+id$ in $d_{xy}$ phase in the experiments of $\kappa$-(BEDT-TTF)${}_{2}\mathrm{X}$.

I thank H. Fukazawa, J. Goryo, H. Kurasawa, H. Nakada, T. Ohama, and Y. Ohta for useful discussions. Parts of numerical calculations were done using the computer facilities of the IMIT at Chiba University, ISSP, and Yukawa Institute.