Hexadecapole fluctuation mechanism for s-wave heavy fermion superconductor CeCu2Si2: Interplay between intra- and inter-orbital Cooper pairs

TL;DR

This paper proposes a new mechanism for s-wave superconductivity in CeCu2Si2, showing that hexadecapole fluctuations driven by vertex corrections induce strong on-site Coulomb attraction, leading to fully-gapped pairing.

Contribution

It introduces a multipole fluctuation theory highlighting hexadecapole fluctuations as the origin of pairing in heavy-fermion superconductors, without phonon involvement.

Findings

Hexadecapole fluctuations mediate attractive pairing.

Superconductivity driven by Coulomb repulsion alone.

Fully-gapped s-wave state explained by multipole interactions.

Abstract

In heavy-fermion superconductors, it is widely believed that the superconducting gap function has sign-reversal due to the strong electron correlation. However, recently discovered fully-gapped s-wave superconductivity in CeCu2Si2 has clarified that strong attractive pairing interaction can appear even in heavy-fermion systems. To understand the origin of attractive force, we develop the multipole fluctuation theory by focusing on the inter-multipole many-body interaction called the vertex corrections. By analyzing the periodic Anderson model for CeCu2Si2, we find that hexadecapole fluctuations mediate strong attractive pairing interaction. Therefore, fully-gapped s-wave superconductivity is driven by pure on-site Coulomb repulsion, without introducing electron-phonon interactions. The present theory of superconductivity will be useful to understand rich variety of the superconducting states in heavy fermion systems.

Tables

Table 1. Table 1: Simple expressions of active EM operators in the Γ 7 ( 1 ) superscript subscript Γ 7 1 \Gamma_{7}^{(1)} - Γ 7 ( 2 ) superscript subscript Γ 7 2 \Gamma_{7}^{(2)} quartet model.

IR ($\Gamma$)	rank (k)	operator ($\hat{Q}$)	matrix
	$0$	$\hat{C}$	${\hat{\sigma }}^0{\hat{\tau }}^0$
$A_1^{+}$	$2$	${\hat{O}}_{20}$	${\hat{\sigma }}^0\left(3{\hat{\tau }}^z+2{\hat{\tau }}^0\right)$
	$4$	${\hat{H}}_0$	${\hat{\sigma }}^0\left(-2.2{\hat{\tau }}^x+2{\hat{\tau }}^z-{\hat{\tau }}^0\right)$
$A_2^{+}$	$4$	${\hat{H}}_z$	${\hat{\sigma }}^z{\hat{\tau }}^y$
$E^{+}$	$2$	${\hat{O}}_{yz},{\hat{O}}_{zx}$	${\hat{\sigma }}^x{\hat{\tau }}^y$, ${\hat{\sigma }}^y{\hat{\tau }}^y$

Table 2. Table 2: Normalized Coulomb interaction for Q 𝑄 Q -multipole channel, U 0 Q subscript superscript 𝑈 𝑄 0 U^{Q}_{0} .

Q	C	$O_{20}$	$H_0$	$H_z$	$O_{yz\left(zx\right)}$
$U_0^Q$	-1.3	-0.18	0.17	0.34	0.27

Equations

$${\hat{H}}_0=\sum _{\mathbf{k}\sigma }{ϵ}_{\mathbf{k}}{c}_{\mathbf{k}\sigma }^{†}{c}_{\mathbf{k}\sigma }+\sum _{\mathbf{k}l\sigma }{E}_{\mathbf{k}l}{f}_{\mathbf{k}l\sigma }^{†}{f}_{\mathbf{k}l\sigma }+\left(V_{\mathbf{k}l\sigma }^{\ast }{f}_{\mathbf{k}l\sigma }^{†}{c}_{\mathbf{k}\sigma }+h.c\right),$$

$$\hat{\chi }(q)=\hat{\varphi }(q)(\hat{1}-u{\hat{U}}^0\hat{\varphi }(q){)}^{-1},$$

$$\hat{\chi }(\mathbf{q},0){\vec{w}}_{\mathbf{q}}^{\Gamma }={\bar{\chi }}_{\mathbf{q}}^{\Gamma }{\vec{w}}_{\mathbf{q}}^{\Gamma },$$

$$\lambda \Delta (k)=\frac{\pi T}{(2\pi {)}^2}\sum _{{ϵ}_m}{\oint }_{\mathrm{FS}}\frac{d\mathbf{p}}{v_{\mathbf{p}}}\frac{\Delta (p)}{|{ϵ}_m|}{V}_{kp}^{\mathrm{sing}},$$

$${\hat{V}}_{kp}^{\chi }=\{u^2{\hat{\Lambda }}_{kp}{\hat{U}}^0\hat{\chi }(k-p){\hat{U}}^0{\hat{\Lambda }}_{kp}^{\prime }{\}}_{\uparrow \uparrow \downarrow \downarrow }-\{{\}}_{\uparrow \downarrow \uparrow \downarrow },$$

$${\hat{V}}_{kp}^U=u\{{\hat{\Lambda }}_{kp}{\hat{U}}^0{\hat{\Lambda }}_{kp}^{\prime }{\}}_{\uparrow \uparrow \downarrow \downarrow }-\{{\}}_{\uparrow \downarrow \uparrow \downarrow },$$

$$\hat{\chi }(q)=\sum _{Q{Q}^{\prime }}{\hat{\chi }}^{Q{Q}^{\prime }}(q),{\hat{\chi }}^{Q{Q}^{\prime }}(q)\equiv {\chi }_q^{Q{Q}^{\prime }}\vec{Q}({\vec{Q}}^{\prime }{)}^{†},$$

$${\hat{V}}_{kp}^{\chi ,Q{Q}^{\prime }}=\{u^2{\hat{\Lambda }}_{kp}{\hat{U}}^0{\hat{\chi }}^{Q{Q}^{\prime }}(k-p){\hat{U}}^0{\hat{\Lambda }}_{kp}^{\prime }{\}}_{\uparrow \uparrow \downarrow \downarrow }-\{{\}}_{\uparrow \downarrow \uparrow \downarrow }.$$

$$\lambda \hat{\Delta }\approx T\sum _k\bar{g}{}^t{\hat{P}}_{\mu \nu }{\hat{G}}^f(k)\hat{\Delta }{\hat{G}}^f(-k){\hat{P}}_{\mu \nu }\theta ({\omega }_c-|{ϵ}_n|).$$

$$G_{lm}^f(k)=G_l^0(k){\delta }_{l,m}+(-1{)}^{l-m}{r}_{\mathbf{k}}{G}_l^0(k)G_m^0(k)G^c(k),$$

$$\hat{G}\hat{\Delta }\hat{G}$$

$$\lambda =g(1+a{)}^2\mathrm{for}(\mu \nu )=(00),(0x),$$

$$\lambda =g(1-a{)}^2\mathrm{for}(\mu \nu )=(0z),(xy),(yy),(zy),$$

$${\hat{H}}_{ff}=\sum _{\mathbf{k}l\sigma }{E}_{\mathbf{k}l}{f}_{\mathbf{k}l\sigma }^{†}{f}_{\mathbf{k}l\sigma }.$$

$$U_0^{Q,Q^{\prime }}=(\vec{Q}{)}^{†}{\hat{U}}^0{\vec{Q}}^{\prime }.$$

$$A_1^{+}$$

$$A_2^{+}$$

$$E^{+}$$

$$A_1^{-}$$

$$A_2^{-}$$

$$E^{-}$$

$$X_{\alpha \beta }^{AL1}(q)=\frac{T}{2}\sum _{{\alpha }^{\prime }{\alpha }^{\prime \prime }{\beta }^{\prime }{\beta }^{\prime \prime }p}{C}_{{\alpha }^{\prime }{\beta }^{\prime \prime }}^{\alpha }(q,p)V_{{\alpha }^{\prime }{\beta }^{\prime }}(p-q)$$

$$\times V_{{\alpha }^{\prime \prime }{\beta }^{\prime \prime }}(p)C_{{\beta }^{\prime }{\alpha }^{\prime \prime }}^{\beta \ast }(\bar{q},\bar{p}),$$

$$C_{ABCD}^{EF}(q,p)\equiv -T\sum _k{G}_{AF}^f(k-q)G_{EC}^f(k)G_{DB}^f(k-p),$$

$$X_{\alpha \beta }^{AL2}(q)$$

$$C_{ABCD}^{{}^{\prime }EF}(q,p)\equiv -T\sum _k{G}_{BF}^f(k-q)G_{ED}^f(k)G_{CA}^f(k-q+p),$$

$${\tilde{C}}_{ABCD}^{{}^{\prime }EF}(q,p)\equiv -T\sum _k{G}_{AE}^f(k+q)G_{FC}^f(k)G_{DB}^f(k+q-p).$$

$$({\hat{\Lambda }}_{k{k}^{\prime }}{)}_{L{L}^{\prime }M{M}^{\prime }}={\delta }_{LM}{\delta }_{L^{\prime }{M}^{\prime }}+({\hat{L}}_{k{k}^{\prime }}{)}_{L{L}^{\prime }M{M}^{\prime }}.$$

$$({\hat{L}}_{k{k}^{\prime }}{)}_{L{L}^{\prime }M{M}^{\prime }}$$

$$B_{ABCDEF}^{M{M}^{\prime }}(q,p,k^{\prime })=G_{AB}^f(k^{\prime }-p)$$

$$\times \left\{C_{CDEF}^{{}^{\prime \prime }M{M}^{\prime }}(q,p)+C_{EFCD}^{{}^{\prime \prime }M{M}^{\prime }}(q,q+p)\right\}$$

$$C_{CDEF}^{{}^{\prime \prime }AB}(q,p)\equiv -T\sum _{k^{\prime }}{G}_{CA}^f(k^{\prime }+q)G_{BF}^f(k^{\prime })G_{ED}^f(k^{\prime }-p).$$

Full text

Heavy fermion (HF) systems exhibit wide variety of unconventional superconductivities Stewart ; Peid ; Maple . For example, antiferro- and ferro- magnetic dipole (rank 1) fluctuations mediate interesting pairing states, such as $d$-wave singlet pairing in Ce$M$In5 ($M$=Rh,Co,Ir) Izawa-115 and triplet pairing in UCoGe Ishida-UCoGe . Since the magnetic dipole fluctuations mediate repulsive pairing interaction, the superconducting gap function inevitably has sign-reversal Moriya ; Yamada ; Kontani-rev ; Scalapino ; Takimoto-SC . However, there are many pairing states in HF systems that cannot be understood based by the rank 1 fluctuations mechanism. In HF systems, it is noteworthy that higher-rank ($r\geq 2$) multipole operators are also active thanks to the strong spin-orbit interaction (SOI), and therefore rich multipole physics emerges. Although higher-rank multipole fluctuations in principle cause exotic pairing states, theoretical studies have not been performed enough.

CeCu2Si2 is a famous HF superconductor near the magnetic criticality Ste-122 ; Yuan-122 ; Ishida ; Hol-122 , and recently reported fully-gapped structure in CeCu2Si2 attract considerable attention Kit1-122 ; Kit2-122 ; Yama-122 ; Steglich . The absence of nodes is confirmed by the measurements of the specific heat, thermal conductivity and penetration depth for $T\ll T_{\rm c}$. In addition, robustness of $T_{\rm c}$ against randomness strongly indicates the plain $s$-wave state without any sign-reversal Yama-122 . Theoretically, magnetic multipole (MM) ($r=1,3,5$) fluctuations will cause sign-reversing pairing states Ikeda-122 . Therefore, electric multipole (EM) ($r=2,4$) fluctuations that give attractive pairing interaction would be important in CeCu2Si2, whereas the microscopic origin of EM fluctuations is unknown.

The minimum theoretical model of CeCu2Si2 is the four-orbital ($J_{z}=\pm 5/2,\pm 3/2$) periodic Anderson model (PAM) with on-site Coulomb interaction. However, if we apply the random-phase-approximation (RPA) to this model, none of EM fluctuations develop. This negative result indicates the significance of the vertex corrections (VCs), which represent the many-body effects beyond the RPA. Recently, it was revealed that higher-rank multipole fluctuations develop cooperatively due to the Aslamazov-Larkin (AL) type $\chi$-VC, which is the VC for the susceptibility, in the study of multipole order in CeB6 Tazai-CeB6 . Physically, the AL-VC gives strong interference between EM and MM fluctuations. Also, the attractive pairing interaction (such as phonon-mediated interaction) is strongly magnified by the $U$-VC, which is the VC for the electron-boson coupling in the gap equation Tazai-PRB2018 . Considering these VCs properly, mysterious plain $s$-wave superconductivity in CeCu2Si2 may be understood in terms of the EM fluctuation mechanism, even if the $e$-ph interaction is absent.

In this paper, we develop a theory of multipole fluctuation mediated superconductivity in HF systems based on the multiorbital PAM. Due to the AL-VC for susceptibility ($\chi$-VC), strong quadrupole and hexadecapole fluctuations develop even in the absence of $e$-ph interaction. In CeCu2Si2, the hexadecapole fluctuations mediate strong attractive pairing interaction, and it is magnified by the AL-VC for the electron-boson coupling ($U$-VC) in the gap equation. Thus, fully-gapped $s$-wave state is caused by the hexadecapole fluctuations against strong on-site repulsive Coulomb interaction. The present pairing mechanism may be significant to understand various HF superconductors.

Now, we introduce a two-dimensional $J=5/2$ PAM for CeCu2Si2. According to the LDA+DMFT study LDADMFT_multiporbital , the following $f$-electron states in $J_{z}$-basis are important near the Fermi level: $|f_{1},\Sigma\rangle=|\mp\frac{5}{2}\rangle$ and $|f_{2},\Sigma\rangle=-|\pm\frac{3}{2}\rangle$, where $\Sigma=\pm$ denotes pseudo-spin Tazai-PRB2018 ; RIXS . The kinetic term of the $\Gamma_{7}^{(1)}$-$\Gamma_{7}^{(2)}$ quartet PAM is given by

[TABLE]

where $c^{\dagger}_{{{\bm{k}}}\sigma}$ ($f^{\dagger}_{{{\bm{k}}}l\sigma}$) is a creation operator for $s\ (f_{l})$-electron with momentum ${{\bm{k}}}$. Here, we put $\Sigma=\sigma$ since the pseudo-spin is conserved in the present PAM Tazai-PRB2018 . We set $\epsilon_{{{\bm{k}}}}=2t_{ss}(\cos k_{x}+\cos k_{y})+\epsilon_{0}$ and $E_{{{\bm{k}}}l}=E^{f}_{l}-(-1)^{l}\delta E_{{{\bm{k}}}}$ ($l=1,2$). Here, $\delta E_{{{\bm{k}}}}$ is given by small $f$-$f$ hopping integrals ($|\delta E_{{{\bm{k}}}}|<0.12|t_{ss}|$) as we explain in the supplemental material (SM) ASM . $V_{{{\bm{k}}}l\sigma}$ is the $f$-$s$ hybridization term between the nearest sites, given as $V_{{{\bm{k}}}l\sigma}=(-1)^{l}t_{sf}^{l}(\sin k_{y}-i\sigma\sin k_{x})$ Tazai-PRB2018 . To make the analysis simple, we set $E_{1}^{f}=E_{2}^{f}\equiv E^{f}$ and $t_{sf}^{1}=t_{sf}^{2}\equiv t_{sf}$. Then, the relation $D_{1}({\epsilon})=D_{2}({\epsilon})$ holds, where $D_{l}({\epsilon})$ is the density of states (DOS) of $f_{l}$-electrons. This is consistent with the relation $D_{1}(0)\approx D_{2}(0)$ given by LDA+DMFT study of CeCu2Si2 LDADMFT_multiporbital ; comment . In the following numerical study, we set $t_{ss}=-1.0$, $E^{f}=0.1$, $\epsilon_{0}=3.0$, $t_{sf}=0.62$, temperature $T=0.045$ and the chemical potential $\mu=-0.143$. Then, $f(s)$-electron number is $n_{f}=0.9$ ($n_{s}=0.3$).

In Fig.1 (a), we show the band structure. $\epsilon=0$ corresponds to the Fermi level. The total band width is $W_{D}\sim 10$ (in unit $|t_{ss}|=1$), and $W_{D}\sim 10$eV in CeCu2Si2 Ikeda-122 . The width of quasi-particle band (=the lowest band) is $W_{D}^{qp}\sim 1$. The Fermi surface (FS) is shown in Fig.1 (b). The anisotropy of $f_{l}$-orbital weight on the FS is introduced by $\delta E_{{{\bm{k}}}l}$, which exists in real HF compounds. We call the present PAM with orbital anisotropy the model A. We will discuss later that the orbital anisotropy is favorable for the $s$-wave state.

We introduce the interaction term $\hat{H}_{U}=u\hat{H}^{0}_{U}$. Here, $\hat{H}^{0}_{U}=\frac{1}{4}\sum_{LL^{\prime}MM^{\prime}}U^{0}_{L,L^{\prime};M,M^{\prime}}f_{L}^{\dagger}f_{L^{\prime}}f_{M}f_{M^{\prime}}^{\dagger}$, where $L=(l,\sigma)$ and $M=(m,\rho)$. $\hat{U}^{0}$ is the $16\times 16$ normalized Coulomb interaction, of which the maximum element is unity Tazai-PRB2018 . The pseudo-spin is conserved in $\hat{H}_{U}$.

The present model belongs to $D_{4h}$ point group. The active irreducible representation (IR) are $\Gamma^{+}=A_{1}^{+},A_{2}^{+},E^{+}$ and $\Gamma^{-}=A_{1}^{-},A_{2}^{-},E^{-}$ Tazai-PRB2018 . In TABLE 1, we show the active EM operators and their approximate pseudo-spin representations. The $4\times 4$ matrix form of each multipole operator ${\hat{Q}}$ is shown in the SM B SM .

From now on, we calculate the $f$-electron susceptibilities. The bare irreducible susceptibility is $\chi_{{\alpha},{\beta}}^{0}(q)=-T\sum_{k}G^{f}_{LM}(k+q)G^{f}_{M^{\prime}L^{\prime}}(k)$, where $q\equiv({{\bm{q}}},\omega_{j})=({{\bm{q}}},2j\pi T)$, ${\alpha}\equiv(L,L^{\prime})$ and ${\beta}\equiv(M,M^{\prime})$. $\hat{G}^{f}$ is the $f$-electron Green function without self-energy Tazai-PRB2018 . To go beyond the RPA, we calculate the AL term for $\chi$-VC, $\hat{X}^{\rm{AL}}$. Its diagrammatic expression and analytic one are respectively given in Fig.2 (a) and in the SM C SM . Since $\chi$-VC is important only for EM susceptibilities, we project out the magnetic channel contribution of $\chi$-VC Yamakawa-FeSe ; Tazai-CeB6 ; Onari-SCVC ; rina2 ; rina1 ; Tsuchiizu . We also drop the MT-type VC since its contribution is small Yamakawa-FeSe ; Tazai-CeB6 ; Onari-SCVC ; rina2 ; rina1 ; Tsuchiizu . Then, the $f$-electron susceptibility in the $16\times 16$ matrix form is given as

[TABLE]

where $\hat{\phi}(q)=\hat{\chi}^{0}(q)+\hat{X}^{\rm{AL}}(q)$ is irreducible susceptibility. To derive the multipole susceptibility, we solve the following eigenvalue equation

[TABLE]

where $\vec{w}^{\Gamma}$ is the eigenvector that belongs to the IR $\Gamma$. It is expressed as $\vec{w}^{\Gamma}=\sum_{Q\in\Gamma}b_{Q}{\vec{Q}}$, where $b_{Q}$ is real coefficient and $\vec{Q}$ is $16\times 1$ vector defined as $(\vec{Q})_{\alpha}\equiv(\hat{Q})_{L,M}$ with ${\alpha}=(L,M)$. Then, the largest eigenvalue $\bar{\chi}^{\Gamma}_{{\bm{q}}}$ gives the multipole susceptibility for the IR $\Gamma$.

In Fig.2 (b), we show the obtained $\bar{\chi}^{\Gamma}_{\rm max}\equiv\max_{{\bm{q}}}\{\bar{\chi}^{\Gamma}_{{\bm{q}}}\}$ for each $\Gamma$. With increasing $u$, all the EM fluctuations strongly develop thanks to the AL-VC. Thus, large EM susceptibilities originate from the interference of MM fluctuations, as discussed in the study of multipole order in CeB6 Tazai-CeB6 . For the EM susceptibilities, the maximum position of $\bar{\chi}^{\Gamma}_{{\bm{q}}}$ for $\Gamma=A_{1}^{+}$ is ${{\bm{q}}}\approx(\pi,\pi)$, whereas that for $\Gamma=A_{2}^{+},E^{+}$ is ${{\bm{q}}}\approx(0,0)$. For the MM susceptibilities, the maximum position for $\Gamma=A_{2}^{-},E^{-}$ is ${{\bm{q}}}\approx(\pi/2,\pi/2)$.

In the next stage, we solve the linearized gap equation with $U$-VC introduced in Ref. Tazai-PRB2018 :

[TABLE]

where $k=({{\bm{k}}},{\epsilon}_{n})=({{\bm{k}}},(2n+1)\pi T)$ and $p=({{\bm{p}}},{\epsilon}_{m})=({{\bm{p}}},(2m+1)\pi T)$. $\Delta(k)$ is the gap function on the FS, $\lambda$ is the eigenvalue, and $v_{{{\bm{p}}}}$ is the Fermi velocity. ${V}^{\rm{sing}}_{kp}$ is the spin singlet paring interaction in band basis, given by the unitary transformation of $\hat{V}_{kp}^{\rm{sing}}\equiv\hat{V}^{\chi}_{kp}+\hat{V}^{U}_{kp}$. Here,

[TABLE]

where $V^{\chi}$ ($V^{U}$) gives the pairing interaction due to fluctuations (Coulomb repulsion). $\hat{\Lambda}_{kp}$ and $(\hat{\Lambda}^{\prime}_{kp})_{LL^{\prime}MM^{\prime}}\equiv(\hat{\Lambda}_{kp})_{M^{\prime}ML^{\prime}L}$ are AL-type $U$-VCs Tazai-PRB2018 . The expression of $\hat{\Lambda}_{kp}$ is given in Ref. Tazai-PRB2018 and SM C SM .

The gap equation is schematically shown in the inset of Fig.3, where the black triangle is the $U$-VC. As explained in Ref. Tazai-PRB2018 , $|\hat{\Lambda}|^{2}\gg 1$ for the electric channel in the presence of MM fluctuations. Therefore, the pairing interaction due to hexadecapole or quadrupole fluctuations is strongly enlarged by $|\hat{\Lambda}|^{2}\gg 1$. As shown in Fig.3, when $u>0.55$, the $d_{x^{2}-y^{2}}$-wave state is replaced with the $s$-wave state mediated by the strong EM fluctuations in Fig. 2 (b). The obtained $s$-wave state is fully gapped without sign reversal, consistently with experiments in CeCu2Si2.

Now, we discuss the origin of the $s$-wave superconductivity. For this purpose, we decompose the susceptibility into the summation of $(Q,Q^{\prime})$-channel as Tazai-PRB2018 ; Tazai-CeB6

[TABLE]

where $\chi^{QQ^{\prime}}_{q}$ is s scalar multipole susceptibility. Then, $\{\vec{Q}\}$ forms a complete non-orthogonal basis: $(\vec{Q})^{\dagger}\vec{Q}^{\prime}$ is unity for $Q=Q^{\prime}$, whereas it is zero when $Q$ and $Q^{\prime}$ belong to different IR. Note that $\bar{\chi}^{\Gamma}_{q}$ is the maximum eigenvalue of the Hermite matrix composed of $\chi^{QQ^{\prime}}_{q}$ with $Q,Q^{\prime}\in\Gamma$.

Then, the $(Q,Q^{\prime})$-fluctuation-induced paring interaction in the band basis, $V^{{\chi},{QQ^{\prime}}}_{kp}$, is given by the unitary formation of

[TABLE]

In Fig.4 (a), we show the EM-fluctuation-mediated interaction averaged on the FS, $V^{\chi,{QQ^{\prime}}}\equiv\int_{\rm FS}dkdpV^{\chi,{QQ^{\prime}}}_{kp}/\int_{\rm FS}dkdp$, for $Q=Q^{\prime}$, together with the total EM-fluctuation interaction $V^{\chi,{\rm EM}}\equiv\sum_{QQ^{\prime}}^{EM}V^{\chi,{QQ^{\prime}}}$. In the present model with $\delta E_{{\bm{k}}}$ (model A), the contribution from the hexadecapole ($H_{0}$) fluctuations in the $A_{1}^{+}$ representation is the largest, while other EM fluctuations are also important. For comparison, we analyze the orbital isotropic model with $\delta E_{{\bm{k}}}=0$, which we call the model I. Surprisingly, in the model I, multipole fluctuations other than $H_{0}$ do not contribute to the $s$-wave pairing, irrespective that all EM ($E^{+},A_{2}^{+},A_{1}^{+}$) susceptibilities develop similarly to Fig. 2 (b) for model A. The FS and its orbital character in each model are shown in Figs.4 (c) and (d). In model I, the orbital weight is perfectly isotropic, whereas the shape of FS is almost model-independent.

Figure 4 (e) shows the $s$-wave pairing interactions $V^{\chi}$ and $V^{U}$ averaged over the FS. Both $V^{\chi}$ and $|V^{U}|$ increase with $u$ in both models due to large $U$-VC for the EM channel Tazai-PRB2018 . In model I, the total pairing interaction $V^{\rm sing}=V^{\chi}+V^{U}$ is always negative (=repulsive), so the $d$-wave state appears. In model A, in contrast, $V^{\rm sing}$ becomes positive with $u$ since not only $H_{0}$ fluctuations, but also other EM fluctuations contribute to the attractive pairing when $\delta E_{{\bm{k}}}\neq 0$. Therefore, the fully-gapped $s$-wave state is realized in model A. As shown in Fig. 3, the eigenvalue $\lambda$ for the $s$-wave state is very large because of the retardation effect as we explain in the SM D SM . In fact, $V^{\rm sing}$ due to EM fluctuations is attractive only for lower frequencies.

Finally, we discuss why all EM fluctuations contribute to the $s$-wave state in model A ($\delta E_{{\bm{k}}}\neq 0$). Since the relation $D_{1}({\epsilon})\approx D_{2}({\epsilon})$ holds even if $\delta E_{{\bm{k}}}$ exists, the obtained EM- and MM-fluctuations are similar in both model A and I. On the other hand, the “inter-orbital pairing $\langle f_{{{\bm{k}}}1\uparrow}f_{-{{\bm{k}}}2\downarrow}\rangle$” is suppressed in model A due to the ${{\bm{k}}}$-dependence of the orbital character on the FS. The absence of inter-orbital pairing is favorable for $s$-wave state as we will discuss later.

One may expect that any EM fluctuations causes the attractive pairing interaction. However, some elements of the EM susceptibility ${\hat{\chi}}^{QQ^{\prime}}({{\bm{q}}},0)=\chi^{QQ^{\prime}}_{q}\vec{Q}(\vec{Q^{\prime}})^{\dagger}$ are negative except for $Q=Q^{\prime}=C$, so the cancellation of pairing interaction may occur. (For example, $\{{\hat{\chi}}^{QQ}({{\bm{q}}})\}_{1\uparrow,2\downarrow;1\uparrow,2\downarrow}<0$ for $Q=O_{yz,zx}$.)

Now, we consider the gap equation when the pairing interaction is given as $\hat{V}^{(\mu\nu)}=\bar{g}\vec{P}_{\mu\nu}(\vec{P}_{\mu\nu})^{\dagger}$, where $\bar{g}>0$ and ${\hat{P}}_{\mu\nu}\equiv\hat{{\sigma}}_{\mu}\hat{\tau}_{\nu}$ ($\mu,\nu=0,x,y,z$). All the EM operators are given by (linear combination of) $\hat{P}_{\mu\nu}$ with $(\mu\nu)=(00),(0x),(0z),(xy),(yy),(zy)$. The gap equation with BCS type cut-off ${\omega}_{c}$ in the orbital basis is

[TABLE]

As explained in Ref. Tazai-PRB2018 , $\hat{G}^{f}$ is expressed as

[TABLE]

where $G^{c}(k)=(i{\epsilon}_{n}-{\epsilon}_{{\bm{k}}}-\sum_{l}r_{{{\bm{k}}}}G_{l}^{0}(k))^{-1}$ is the $c$-Green function, $r_{{{\bm{k}}}}\equiv|V_{{{\bm{k}}}1\uparrow}|^{2}$, and $G_{l}^{0}(k)=(i{\epsilon}_{n}-E_{{{\bm{k}}}l})^{-1}$ is the unhybridized $f$-Green function. We neglect the first term in Eq. (8) since it does not give $-{\rm ln}T$ term in gap equation. Therefore, in model I ($\delta E_{{\bm{k}}}=0$), the relation $\hat{G}\propto\hat{{\sigma}}_{0}(\hat{\tau}_{0}-\hat{\tau}_{x})$ holds since $G_{1}^{0}=G_{2}^{0}$. In model A with large $|\delta E_{{\bm{k}}}|\ (\gg|V_{{{\bm{k}}}1\uparrow}|^{2})$, the relation $\hat{G}\propto\hat{{\sigma}}_{0}\hat{\tau}_{0}$ holds approximately.

Here, we set $\hat{G}\propto\hat{\sigma}_{0}(\hat{\tau}_{0}-a\hat{\tau}_{x})$ and $\hat{\Delta}\propto i\hat{\sigma}_{y}(\Delta_{0}\hat{\tau}_{0}+\Delta_{x}\hat{\tau}_{x})$: $a=1$ $(a=0)$ corresponds to model I (model A with $|\delta E_{{\bm{k}}}|\gg|V_{{{\bm{k}}}1\uparrow}|^{2}$). In this case,

[TABLE]

Then, the eigenvalue of the gap equation is

[TABLE]

where $g=\bar{g}D_{1}(0){\rm ln}({\omega}_{c}/T)$.

In Fig.5, we summarize the eigenvalue $\lambda$ for each EM pairing interaction, in the case of $a=0$ (intra orbital Cooper pair) and $a=1$ (intra+inter orbital Cooper pair). We note that $\hat{P}_{0z}\propto\hat{O}_{20}-2\hat{C}$ and $\hat{P}_{0x}\propto-3\hat{H}_{0}+2\hat{O}_{20}+\hat{C}$. In case of $a=0$, all EM fluctuations contribute to the pairing. In case of $a=1$, however, only $\hat{P}_{0x}$ and $C$ channels contribute to the pairing. In the present PAM, charge ($C$) fluctuations are small, so they do not contribute to the pairing. Since $\hat{P}_{0x}$ is included only in $H_{0}$ hexadecapole, the $H_{0}$ fluctuations give dominant $s$-wave pairing interaction. To summarize, the pairing interaction increases if the inter orbital Cooper pairs are killed by finite $|\delta E_{{\bm{k}}}|$, so the numerical results in Fig.4 are well understood.

In summary, we studied the multipole fluctuation mediated superconductivity in HF systems based on the $\Gamma_{7}^{(1)}$-$\Gamma_{7}^{(2)}$ quartet PAM. Due to the AL-type $\chi$-VC, strong quadrupole and hexadecapole fluctuations develop, and the resultant attractive interaction is enlarged by the $U$-VC in the gap equation. In CeCu2Si2, $H_{0}$ hexadecapole fluctuations mediate strong attractive pairing interaction. The $s$-wave state is further stabilized by introducing small $\delta E_{{{\bm{k}}}l}$, by which the inter-orbital Cooper pairs are killed. Moreover, if we introduce the $e$-ph interaction, both $U$-VCs ($\chi$-VC and $U$-VC) and $e$-ph interaction would enlarge $s$-wave $T_{\rm c}$ cooperatively Tazai-PRB2018 ; Razaf ; Ohkawa ; Nagaoka ; Kontani-RPA . The present pairing mechanism may be significant to understand various HF superconductors.

There are many important future issues, such as the self-energy effect Haule ; LDADMFT_multiporbital ; KotGeo-DMFT ; KotVol-DMFT ; Held ; DMFT-HF and verificatoin of the multiorbital nature of the FS comment ; Fulde-comment ; LDADMFT_multiporbital ; Hat-122 . Also, $P$-induced second superconducting phase of CeCu2Si2 is an important issue Hol-122 .

.1 A: Model Hamiltonian

In this section, we introduce the $f$-$f$ hopping integrals, by which the $f_{l}$-orbital weight has momentum-dependence on the FS. The obtained $f$-electron kinetic term is S-Tazai-PRB2018

[TABLE]

Here, we set $E_{{{\bm{k}}}1}\equiv E_{1}+\delta E_{{{\bm{k}}}}$ and $E_{{{\bm{k}}}2}\equiv E_{2}-\delta E_{{{\bm{k}}}}$. To reproduce the ${{\bm{k}}}$-dependent $\delta E_{{{\bm{k}}}}$ shown in Fig.S1, we introduce from the first to fifth neighbor hopping integrals according to Refs.S-Yamakawa-FeSe ; S-Tazai-PRB2018 . The obtained momentum-dependence of $f_{l}$-orbital weight on the FS is shown in Fig. 1 (b) in the main text.

As discussed in Ref. S-Tazai-PRB2018 , the RPA susceptibility is insensitive to $\delta E_{{{\bm{k}}}}$ since the $f_{l}$-orbital DOS, $D_{l}({\epsilon})$, is independent of $\delta E_{{{\bm{k}}}}$. In the present study, we verified that both $\chi$-VC and $U$-VC are also insensitive to $\delta E_{{{\bm{k}}}}$.

In HF systems, the quadrupole susceptibility remains small within the RPA. To understand this result, we examine the $(Q,Q^{\prime})$ component of normalized Coulomb interaction:

[TABLE]

TABLE 2 shows the diagonal component $U^{Q}_{0}\equiv U^{Q,Q}_{0}$. Since $U^{Q}_{0}$ for the EM channels is much smaller than that for the MM channels, the EM susceptibilities are small within the RPA. Nonetheless of this fact, EM susceptibilities strongly develop by considering the AL-VC, since $X^{\rm AL}$ for the EM channel becomes large when moderate MM fluctuations exist.

.2 B: Pseudospin representation of multipole operators

Here, we list the multipole operators $\hat{Q}$ in the present CeCu2Si2 model, which were already explained in Ref. S-Tazai-PRB2018 . The EM (even-rank) operators in the $4\times 4$ matrix form are expressed as

[TABLE]

The MM (odd-rank) operators are given by

[TABLE]

where $\hat{\sigma}^{\mu}$ and $\hat{\tau}^{\mu}$($\mu=x,y,z$) are Pauli matrices for the pseudo-spin and orbital basis, respectively. $\hat{\sigma}^{0}$ and $\hat{\tau}^{0}$ are identity matrices.

The row and column of the Hermite matrix $\hat{Q}$ for each operator is given as $L=(l,{\sigma})$, where $l=1,2$ represents the $f$-orbital and ${\sigma}=\uparrow,\downarrow$ represents the pseudo spin. In the main text, we also introduce the vector representation defined as $(\vec{Q})_{\alpha}=(\hat{Q})_{L,L^{\prime}}$, where ${\alpha}=(L,L^{\prime})$.

.3 C: Analytic expressions of vertex corrections

From now on, we introduce the analytic expressions of $\chi$-VC S-Tazai-CeB6 and $U$-VC S-Tazai-PRB2018 due to AL diagrams. First, we discuss the $\chi$-VCs, whose diagrammatic expressions are shown in Fig. 2 (a) in the main text. The expression for the AL1 term is given as

[TABLE]

where $p\equiv({{\bm{p}}},\omega_{j})$, $\bar{p}\equiv({{\bm{p}}},-\omega_{j})$, and $\hat{V}(q)\equiv u^{2}\hat{U}^{0}\hat{\chi}(q)\hat{U}^{0}+u\hat{U}^{0}$ is the dressed interaction given by the RPA. The three-point vertex in Eq. (S9) is given as

[TABLE]

where $\hat{G}^{f}$ is the $f$-electron Green function. Also, the expression for the AL2 term is given as

[TABLE]

where

[TABLE]

The total $\chi$-VC is given by $\hat{X}^{\rm AL}=\hat{X}^{AL1}+\hat{X}^{AL2}$, by subtracting the double counting second order diagrams of order $u^{2}$.

Next, we explain the $U$-VC in the gap equation. It is given as

[TABLE]

In the main text, we calculate the AL diagrams for ${\hat{L}}_{kk^{\prime}}$. It is expressed as

[TABLE]

where

[TABLE]

and

[TABLE]

.4 D: Gap equation and retardation effect

Here, we comment on the retardation effects. In Fig.S2, we show the obtained paring interaction on the FS defined as $V^{\rm{sing}}_{\rm max}(\omega_{j})\equiv\max_{{{\bm{k}}},{{\bm{k}}}^{\prime}}\{V^{\rm{sing}}_{({{\bm{k}}},\pi T)({{\bm{k}}}^{\prime},\pi T+\omega_{j})}\}$. The paring interaction is attractive (positive) at $\omega_{j}=0$, whereas it becomes to repulsion for $\omega_{j}>0$. For this reason, the gap function defined as $\Delta(\epsilon_{n})\equiv\max_{{{\bm{k}}}}\{\Delta({{\bm{k}}},\epsilon_{n})\}$ shows the sign-change as the function of ${\epsilon}_{n}$, as shown in the inset of Fig.S2. This is a hallmark of the retardation effects due to the strong $\omega_{j}$-dependence of the EM (even-rank) fluctuation. Since the depairing due to direct Coulomb interaction is reduced by the retardation effect, the fully-gapped $s$-wave superconductivity can be stabilized in HF systems.