Prediction of pseudogap formation due to $d$-wave bond-order in organic superconductor $\kappa$-(BEDT-TTF)$_2$X

Rich hidden unconventional orders with pseudogap formation, such as the inter-site bond-order (BO), attract increasing attention in condensed matter physics. Here, we investigate the hidden order formation in organic unconventional superconductor $\kappa$-(BEDT-TTF)$_2$X. We predict the formation of $d$-wave BO at wavelength $q=Q_B=(\delta,\delta)$ ($\delta=0.38\pi$) for the first time, based on both the functional renormalization group (fRG) and the density-wave equation theories. The origin of the BO is the quantum interference among antiferromagnetic spin fluctuations. This prediction leads to distinct pseudogap-like reduction in the NMR $1/T_1$ relaxation rate and in the density-of-states, consistently with essential experimental reports. The present theory would be applicable for other strongly correlated metals with pseudogap formation.

In this paper, we discuss the occurrence of hidden charge DW orders in κ-(BEDT-TTF) 2 X by using both the functional renormalization group (fRG) and the DW equation methods. We predict the d-wave BO order formation at wavelength q = Q B = (δ, δ) (δ = 0.38π) in κ-(BEDT-TTF) 2 X for the first time. The origin of the BO is novel quantum interference among paramagnons with Q ± S ≈ (π ± δ/2, π ± δ/2). The predicted distinct reductions in the spin fluctuation strength and in the DOS at BO transition temperature naturally explain the NMR and Raman measurements.
To analyze the quantum interference in lowdimensional metals, the fRG method is very powerful and reliable, since various DW instabilities (both particleparticle and particle-hole channels) are treated on the same footing [41,46,47,[55][56][57][58][59]. Using the fRG, unconventional DW states in cuprates and ruthenates have been analyzed satisfactorily [47,58]. The fRG method is suitable to analyze the many-body electronic states, and we can obtain reliable results by making careful comparison with the diagrammatic calculation using the DW equation method [47,[51][52][53].
In the following numerical study, we set the energy unit |t| = 1, and put the temperature T = 0.05 and the electron filling n = 1 (µ = 0.55). The band structure and the Fermi surface (FS) are presented in Figs. 1 (b) and (c), respectively. The patch indices (1 ∼ 64) are shown on the ellipsoid electron pockets. The total band width is W D ∼ 10 (in unit |t| = 1), and |t| corresponds to 0.05eV since W D ∼ 0.5eV experimentally [5,60].
In the present study, we analyze the model by applying the RG+cRPA method developed in Refs. [46,47,58,59], which is an efficient hybrid method between the fRG and the random-phase-approximation (RPA). Here, we introduce the higher-energy cutoff Λ 0 (= 2) and the logarithmic energy mesh Λ l = Λ 0 e −l with l ≥ 0. Then, the effective four-point vertexΓ is derived from the oneloop RG equation in Fig. 1 (d), by including the onshell contribution Λ l+dl < |ǫ k − µ| ≤ Λ l step-by-step. In strongly correlated systems,Γ strongly deviates fromΓ 0 . In this procedure, the Brillouin zone (BZ) is divided into N p patches; see the Supplemental Materials (SM) A [62]. Below, we set Λ 0 = 2 and N p = 64. (Numerical results are not sensitive to the choices of Λ 0 and N p .) We analyze the higher energy region |ǫ k − µ| ≥ Λ 0 , where the vertex corrections (VCs) are less important, based on the constrained RPA (cRPA) with high numerical accuracy using fine k-meshes. By using the RG+cRPA method, we can perform reliable calculations [46,47,58].
Here, we show numerical results for the dimer Hubbard model. First, we calculate various DW susceptibilities at ω = 0 using the RG+cRPA method. The spin (charge) susceptibility with the form factor f x (k) is [47] where A s(c) We calculate the spin and charge susceptibilities with the following form factors cos k y and f xy = 2 sin k x sin k y . As a result, we find that the spin susceptibility χ S (q) ≡ χ s 1 (q) and the d x 2 −y 2 -wave BO susceptibility χ BO (q) ≡ χ c x 2 −y 2 (q) strongly develop. Other susceptibilities remain small in the present study.
In Figs. 2 (a) and (b), we plot q-dependences of χ S (q) and χ BO (q) at U = 3.5. Strong spin fluctuations develop at q = Q S , consistently with the previous RPA and FLEX analyses [5]. In addition, we reveal the development of χ BO (q) at q = Q B ≈ (3π/8, 3π/8) in addition to q = (0, 0). The obtained strong bond fluctuations originate from the VCs that are dropped in the RPA.
The χ BO (q) strongly develops by increasing U . Figure  2 (c) shows the RG flow of the susceptibilities in the case of U = 3.54. In this case, the bond susceptibility exceeds the spin one after completing the renormalization. We see that χ S (Q S ) starts to increase in the early stage of the renormalization, by reflecting the major nesting of the FS at q = Q S . Next, χ BO (Q B ) starts to increase for l 3, and it exceeds χ S (Q S ) at l ∼ 4. Finally, χ BO (0) starts to increase for l 4 (Λ l 0.037), because the renormalization of the Pauli (q = 0) susceptibility occurs only for Λ l T . All susceptibilities saturate for l 8 (Λ l 0.7 × 10 −3 ). The final results in Figs. 2  We also calculate the spin-singlet SC susceptibility [47] (2) where ∆(k) is an even parity gap function, which is uniquely determined so as to maximize χ SC under the constraint 1 N k |∆(k)| 2 δ(ǫ k − µ) = 1 [47]. We show the obtained optimized gap at U = 3.54 in Fig. 2 (d).
RG equation by dropping the contribution from the pp channel in the RG equation for Γ, by introducing an additional cutoff energy only for the pp channel; ω pp c (> ω c ). Here, we set ω pp c = 10T to suppress the SC fluctuations selectively. The obtained RG flows of the susceptibilities at U = 3.82 are shown in Fig. 3 (a). We see that χ S (Q S ) starts to increase in the early stage, due to the ph channels in the RG equations in Fig. 1 (d). Next, χ BO (Q B ) also increases to follow the increment of χ S (Q S ), similarly to Fig. 2 (c). This result strongly indicates that the BO fluctuations are driven by the spin fluctuations. Note that χ BO (Q B ) exceeds χ S (Q S ) by setting U = 3.86 even in the case ω pp c = 10T . Next, we discuss the SC pairing vertex function V SC (p, p ′ ) = 3 2 Γ s (p, p ′ , −p ′ , −p) − 1 2 Γ c (p, p ′ , −p ′ , −p) − Γ 0s given by the RG+cRPA method [59]. Due to large ω pp c (= 10T ), the obtained V SC (p, p ′ ) becomes "irreducible with respect to the pp channel" below ω pp c . Then, V SC (p, p ′ ) gives the "pairing interaction in the SC gap equation" with the BCS cutoff energy ω BCS = ω pp c . The obtained RG flow of V SC (p, p ′ ) is shown in Fig. 3 (b). The large repulsion for p − p ′ ≈ Q S is apparently given by the spin fluctuations. Interestingly, we find that the attraction for p − p ′ ≈ Q B is caused by the BO fluctuations in Fig. 3 (a). The present result indicates that both χ BO (Q B ) and χ S (Q S ) cooperatively work as the pairing glue of the d x 2 −y 2 -wave state, as understood in Fig. 1 (c). The obtained gap structure is very similar to Fig. 2 (d). Therefore, the d-wave BO fluctuations in the single-orbital Hubbard model can mediate large attractive pairing interaction.
Next, we explain that the BO fluctuations originate from the quantum interference between paramagnons, which is described by the Aslamazov-Larkin (AL) quantum process. For this purpose, we analyze the following DW equation [48,[51][52][53]: where λ DW q is the eigenvalue that represents the charge channel DW instability at wavevector q. Here, p ± ≡ p ± q/2, and k ≡ (k, ǫ n ) and p ≡ (p, ǫ m ) (ǫ n , ǫ m are fermion Matsubara frequencies). The eigenfunction f q (k) gives the form factor. The corresponding DW susceptibility is In the present model, the obtained DW state corresponds to the d-wave BO.
The kernel function I c is given by the Ward identity −δΣ/δG, which is composed of one single-paramagnon exchange term and two double-paramagnon exchange ones: The former and the latter are called the Maki-Thompson (MT) term and the AL terms; see Fig. 4 (a). (Each wavy line is proportional to χ S (q).) Among four terms in I c q , the AL terms are significant for α S 1, and they give spin-fluctuation-driven nematic orders in cuprates and Fe-based superconductors [48,49].  reaches almost unity at q = Q B , which is consistent with χ BO (q) given by the RG+cRPA in Fig. 2  (b). The corresponding form factor in Fig. 4 (c) possesses the d x 2 −y 2 -wave symmetry. Therefore, strong d x 2 −y 2wave BO susceptibility at q = Q B by the RG+cRPA method in Fig. 2 (b) is well reproduced by the DW equa-tion. The origin of the strong BO instability at q = Q B is the AL terms in Fig. 4 (a) that represent the quantum interference among paramagnons at Q ± S ≈ (π, π)±Q B /2. (Note that χ S (q) given by the RPA is similar to that the fRG result in Fig. 2 (a).) The paramagnon-interference mechanism can generate both the ferro-BO instability (at q = Q ± S − Q ± S = 0) and the incommensurate-BO one (at q = Q + S − Q − S = Q B ). This mechanism causes the ferro-BO states in both Febased and cuprate superconductors according to the DW equation analysis [48,[51][52][53]. In the present dimer Hubbard model, in contrast, the ferro-BO fluctuations remain small in the DW equation analysis. This is also true in the fRG analysis with ω pp c = 10T shown in Fig. 3. These results indicate that the paramagnon interference mechanism alone is not sufficient to establish large χ BO (0) in Fig. 2 (b). Therefore, we conclude that large χ BO (0) is caused by the spin and SC fluctuations cooperatively, since the AL processes by SC fluctuations can cause the ferro-BO fluctuations according to Ref. [58].
Finally, we discuss the band-folding and hybridization gap due to the BO with q = Q B . Figure 4 (d) shows the Fermi arc structure obtained for f max ≡ max k {f QB (k)} = 0.1. Here, the folded band structure under the BO at q = Q B is "unfolded" into the original Brillouin zone [63] to make a comparison with ARPES experiment. The resultant pseudogap in the DOS is shown in the inset of Fig. 4 (d), which is consistent with the STM study [11]. The BO leads to significant reduction of the spin fluctuation strength, so the obtained 1/T 1 T ∝ q,α,β Im χ s α,β (q, ω)/ω In the fRG study, the parquet VCs are generated by considering all ph, ph' and pp channels in Fig. 1 (d). On the other hand, in the DW equation study, the VCs are limited to MT and AL terms, whereas their frequency dependences are calculated correctly. Both theoretical methods lead to the emergence of the same d-wave bond order shown in Fig. 2 (e).
In summary, we predicted the emergence of the d-wave BO at wavevector q = Q B = (0.38π, 0.38π) in κ-(BEDT-TTF) 2 X, due to the interference between paramagnons with Q ± S ≈ (π, π) ± Q B /2. The BO is derived from both fRG method and the DW equation method. The BO transition leads to distinct pseudogap behaviors in the NMR 1/T 1 relaxation rate and in the DOS, consistently with many experimental reports at T ≈ T * . As we show in the SM B [62], very similar numerical results are obtained in the case of t ′ /t = 0.7. Thus, the d-wave BO will be ubiquitous in κ-(BEDT-TTF) 2 X. The present theory would be applicable for other strongly correlated metals with pseudogap formation.
In the main text, we studied the dimer Hubbard model, which is the simplest effective model for κ-(BEDT-TTF) 2 X, by applying the RG+cRPA method. This is a useful hybrid method between the fRG theory and the RPA developed in Ref. [1]. The RG+cRPA method is applicable for systems with complex valence-band structure, even when the conventional patch RG method is not applicable. Here, we solve the one-loop RG equation for the four-point vertex Γ in Fig. 1 (d) inside the energy region |ǫ k − µ| ≤ Λ 0 , by applying the logarithmic energy mesh Λ l = Λ 0 e −l with l ≥ 0. In the framework of the N p -patch RG, the lower-energy region in the Brillouin zone (BZ) is divided into N p patches, as shown in Fig.  S1 with Λ 0 = 2 and N p = 32. (Note that the numerical study in the main text is done for N p = 64.) By solving the RG equation, we can calculate the vertex corrections (VCs), which are the many-body effects that are dropped in the RPA.  at q = Q B is very similar to the d x 2 −y 2 -wave form factor for t ′ /t = 0.5 obtained in Fig. 4 (c) in the main text. These results are essentially similar to the results by the RG+cRPA method in Fig. S2, except for the absence of the peak at q = 0 in Fig. S3 (b).
In summary, the development of the d-wave BO susceptibility at q = Q B , χ BO (Q B ), is confirmed by both the fRG theory and the DW equation theory in the dimer Hubbard model, in the cases of t ′ /t = 0.5 and 0.7. This result is derived from the AL-type VCs, which are neglected in the RPA. Since λ DW q=0 in Fig. S3 (b) remains small, the enhancement of χ BO (0) in Fig. S3 (c) originates from the spin and SC fluctuations cooperatively, as we discussed in the main text.
Finally, it should be stressed that the ferro-BO can be induced by the paramagnon-interference mechanism in general systems. In fact, strong ferro-BO fluctuations observed in both Fe-based and cuprate superconductors are satisfactorily reproduced by the DW equation analysis [5][6][7], meaning that the ferro-BO fluctuations originate solely from the spin fluctuations.