Multipole fluctuation theory for heavy fermion systems: Application to multipole orders in CeB6

TL;DR

This paper develops a multipole fluctuation theory for heavy fermion systems, revealing how higher-order many-body effects induce inter-multipole couplings that explain complex orders like the antiferro-quadrupole order in CeB6.

Contribution

It introduces a novel theoretical framework analyzing multipole fluctuations via vertex corrections, elucidating the origin of multipole orders in heavy fermion materials.

Findings

Magnetic, quadrupole, and octupole fluctuations develop cooperatively.

Antiferro-quadrupole order is driven by interference between magnetic-multipole fluctuations.

Inter-multipole coupling mechanism may explain hidden orders in heavy fermion systems.

Abstract

In heavy fermion systems, multipole degrees of freedom make possible the emergence of rich phenomena, such as hidden orders and superconductivities. However, many of them remain unsolved since the origin of higher-rank multipole interaction is not well understood. Among these issues, we focus on the quadrupole order in CeB$_6$, which is a famous multipolar heavy fermion system actively studied for decades. We analyze the multiorbital periodic Anderson model for CeB$_6$, and find that both magnetic, quadrupole, and octupole fluctuations develop cooperatively due to the strong inter-multipole coupling given by higher-order many-body effects, called the vertex corrections. It is found that the antiferro-quadrupole order in CeB$_6$ is driven by the interference between magnetic-multipole fluctuations. The discovered inter-multipole coupling mechanism is a potential origin of various hidden orders in various heavy fermion systems.

Tables

Table 1. Table 1: IRs and 16-type active multipole operators of D 4 ​ h subscript 𝐷 4 ℎ D_{4h} point group. Operator with rank k 𝑘 k corresponds to 2 k superscript 2 𝑘 2^{k} -pole.

IR ($\Gamma$)	rank (k)	Operator ($\hat{Q}$)	IR in $H_z$
${\Gamma }_1^{+}$	$0$	$\hat{1}$	${\Gamma }_1$
	$2$	${\hat{O}}_{20}$
${\Gamma }_3^{+}$	$2$	${\hat{O}}_{22}$	${\Gamma }_3$
${\Gamma }_4^{+}$	$2$	${\hat{O}}_{xy}$	${\Gamma }_4$
${\Gamma }_5^{+}$	$2$	${\hat{O}}_{yz},{\hat{O}}_{zx}$	${\Gamma }_5$
${\Gamma }_2^{-}$	$1$	${\hat{J}}_z$	${\Gamma }_1$
	$3$	${\hat{T}}_{z\alpha }$
${\Gamma }_3^{-}$	$3$	${\hat{T}}_{xyz}$	${\Gamma }_4$
${\Gamma }_4^{-}$	$3$	${\hat{T}}_{z\beta }$	${\Gamma }_3$
${\Gamma }_5^{-}$	$1$	${\hat{J}}_x$,${\hat{J}}_y$	${\Gamma }_5$
	$3$	${\hat{T}}_{x\alpha }$,${\hat{T}}_{y\alpha }$
	$3$	${\hat{T}}_{x\beta }$,${\hat{T}}_{y\beta }$

Table 2. Table 2: Normalized Coulomb interaction U 0 Q subscript superscript 𝑈 𝑄 0 U^{Q}_{0} . U 0 Q , Q ′ = 0 subscript superscript 𝑈 𝑄 superscript 𝑄 ′ 0 0 U^{Q,Q^{\prime}}_{0}=0 for Q ≠ Q ′ 𝑄 superscript 𝑄 ′ Q\neq Q^{\prime} except for U 0 J μ , T μ α = 0.58 subscript superscript 𝑈 subscript 𝐽 𝜇 subscript superscript 𝑇 𝛼 𝜇 0 0.58 U^{J_{\mu},T^{{\alpha}}_{\mu}}_{0}=0.58 ( μ = x , y , z 𝜇 𝑥 𝑦 𝑧 \mu=x,y,z ).

Q	1	$O_{20\left(22\right)}$	$O_{xy(yz,zx)}$	$T_{xyz}$	$J_{z(x,y)}$	$T_{z(x,y)}^{\alpha }$	$T_{z(x,y)}^{\beta }$
$U_0^Q$	-2.4	0.50	0.63	0.81	1.03	0.94	0.94

Equations

$$|f_1\Sigma ⟩=\sqrt{\frac{5}{6}}|\mp \frac{5}{2}⟩+\sqrt{\frac{1}{6}}|\pm \frac{3}{2}⟩,$$

$$V_{\mathbf{k}{f}_l\uparrow }$$

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

$$u{\hat{U}}^0\hat{\varphi }(\mathbf{q},0){\vec{w}}^{\Gamma }(\mathbf{q})={\alpha }^{\Gamma }(\mathbf{q}){\vec{w}}^{\Gamma }(\mathbf{q}).$$

$${\chi }^{Q,Q^{\prime }}(q)=(\vec{Q}{)}^{†}\hat{\chi }(q){\vec{Q}}^{\prime }.$$

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

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

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

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

$$X_{O_{xy}}(q)\equiv ({\vec{O}}_{xy}{)}^{†}\hat{X}(q){\vec{O}}_{xy},$$

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

$$X_{O_{xy}}^{AL1,Q{Q}^{\prime }}(q)\equiv \frac{T}{2}\sum _p{v}_p^Q{v}_{p-q}^{Q^{\prime }}{\Lambda }_{q,p}^{O_{xy}Q{Q}^{\prime }}({\Lambda }_{\bar{q},\bar{p}}^{O_{xy}{Q}^{\prime }Q}{)}^{\ast },$$

$${\Lambda }_{q,p}^{O_{xy}Q{Q}^{\prime }}\equiv \sum _{\alpha }({\vec{O}}_{xy}{)}_{\alpha }^{\ast }({\vec{Q}}^{\prime }{)}^{†}{\hat{\Lambda }}^{\alpha }(q,p)\vec{Q}.$$

$${ϵ}_{\mathbf{k}}$$

$$V_{\mathbf{k}{f}_1\uparrow }$$

$$V_{\mathbf{k}{f}_2\uparrow }$$

$${\hat{H}}_U=u{\hat{H}}_U^0,$$

$${\hat{H}}_U^0=\frac{1}{4}\sum _{L{L}^{\prime }M{M}^{\prime }}{U}_{L,L^{\prime };M,M^{\prime }}^0{f}_L^{†}{f}_{L^{\prime }}{f}_M{f}_{M^{\prime }}^{†},$$

$${\Gamma }_1^{+}$$

$${\Gamma }_3^{+}$$

$${\Gamma }_4^{+}$$

$${\Gamma }_5^{+}$$

$${\Gamma }_2^{-}$$

$${\Gamma }_3^{-}$$

$${\Gamma }_4^{-}$$

$${\Gamma }_5^{-}$$

$$\hat{Q}/\sqrt{\sum _{L,M}|Q_{L,M}{|}^2}\to \hat{Q}.$$

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

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

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

$$X_{L{L}^{\prime }M{M}^{\prime }}^{MT}(q)=T^2\sum _{p,k,A\sim D}{G}_{LA}(k+q-p)G_{B{L}^{\prime }}(k-p)$$

$$\times G_{DM}(k+q)G_{M^{\prime }C}(k)V_{DACB}(p).$$

$${\chi }^{\mathrm{mag}}(\mathbf{q},{\omega }_l)$$

Full text

Heavy fermion (HF) systems are very interesting platform of exotic electronic states induced by strong Coulomb interaction and spin-orbit interaction (SOI) on $f$-electrons. Magnetic fluctuations cause interesting quantum critical phenomena and superconductivity Coleman ; Moriya ; Yamada ; Kontani-rev ; Monthoux ; Tremblay ; Scalapino ; Takimoto-SC . In addition, higher-rank multipole operators are also active thanks to the strong SOI of $f$-electrons. For this reason, various interesting multipole order and fluctuations, which are absent in transition metal oxides, emerge in HF systems. One of the most famous example is the multipole-order in CeB6: The antiferro-quadrupole order with ${{\bm{q}}}=(\pi,\pi,\pi)$ occurs at $T_{Q}=3.2$K, and magnetic order appears at $T_{N}=2.4$K Kasuya ; Goto ; Sera ; Inosov-review . In addition, antiferro-octupole order is stabilized under weak magnetic field Shiina1 ; Shiina2 ; Shiina3 ; Shiina4 . Thus, various ranks of multipole orders appear simultaneously in the phase diagram of CeB6. This fact indicates that different multipoles are strongly entangled, which would be universal in HF system.

Up to now, multipole orders in CeB6 has been studied actively based on the localized $f$-multipole models Shiina1 ; Shiina2 ; Shiina3 ; Shiina4 ; Sera2 ; Kusunose ; Hanzawa . However, recent ARPES studies ARPES1 ; ARPES2 revealed that the $f$-electron is itinerant for $T\sim T_{Q}$. The characteristic dynamical magnetic susceptibility in CeB6 measured by neutron inelastic scattering Inosov1 ; Inosov2 is explained in the itinerant picture based on the periodic Anderson model (PAM) Thal . If we apply the random-phase-approximation (RPA) for the PAM, however, quadrupole order cannot be obtained. In fact, only odd-rank (=magnetic) multipole fluctuations develop, whereas even-rank (=electric) multipole ones remain small in the RPA Thal ; Ikeda-U ; Tazai-HF . This fact means that the importance of vertex corrections (VCs), which represent the many-body effects ignored in the RPA. The lowest-order VC with respect to fluctuations, called the Maki-Thompson (MT) type VC, gives the rank-5 multipole order in URu2Si2 Ikeda-U . However, the MT-VC does not magnify the even-rank multipole fluctuations. Thus, microscopic origin of quadrupole order, which frequently appears in various compounds, is still unsolved. CeB6 is a suitable platform to construct a theory of multipole order in HF systems.

In Fe-based and cuprate superconductors, significant roles of the Aslamazov-Larkin (AL) VC, which is the second-order VC with respect to fluctuations, attract considerable attention Onari-SCVC ; Onari-FeSe ; Yamakawa-FeSe . The AL-VC describes various spin-fluctuation-driven nematicities, such as orbital order and bond order, that fail to be explained by the RPA. The significance of the AL-VC near the magnetic criticality is confirmed by different theoretical studies, especially by the functional-renormalization-group (fRG) studies Tsuchiizu1 ; Tsuchiizu2 ; Tsuchiizu3 ; Tazai-RG ; Tsuchiizu4 ; Chubukov-RG1 ; Schmalian . In HF systems, phonon-mediated superconductivity can be stabilized by the AL-process for the electron-boson coupling. This mechanism may be responsible for the fully gapped $s$-wave state in CeCu2Si2 Tazai-HF ; Tazai-JPSJ . These findings indicate that the AL-VC plays essential roles in HF systems

In this paper, we study the mechanism of quadrupole order in CeB6 based on the itinerant $f$-electron picture, by considering the AL-VC for multipole susceptibilities. For this purpose, we introduced an effective PAM for CeB6 with $\Gamma_{8}$ quartet $f$-orbital basis. Both ferro- and antiferro-magnetic and octupole fluctuations are induced by the Fermi surface nesting, consistently with recent neutron experiments. Then, antiferro-quadrupole ($O_{xy}$) order is induced by the interference between different magnetic multipole fluctuations. The present multipole fluctuation theory with introducing AL-VC will be applicable for various HF systems.

In HF systems, the DMFT has been applied successfully LDADMFT_multiporbital ; KotVol-DMFT ; Held ; DMFT-HF ; Otsuki-HF . In the DMFT, the irreducible VC is local. Here, we calculate the ${{\bm{k}}}$-dependence (=nonlocality) of AL and MT diagrams accurately in order to evaluate the interference between multipole fluctuations.

Here, we introduce a two-dimensional PAM as an effective model for CeB6. For $f$-electron states, we consider the $\Gamma_{8}$ quartet in $J=5/2$ space Shiina1 :

[TABLE]

where $\Sigma=\pm$ is the pseudo-spin of $f_{l}$-orbital ($l=1,2$). The kinetic term is given by $\hat{H}_{0}=\sum_{{{\bm{k}}}\sigma}\epsilon_{{{\bm{k}}}}c^{\dagger}_{{{\bm{k}}}\sigma}c_{{{\bm{k}}}\sigma}+\sum_{{{\bm{k}}}l\sigma}E_{f}f^{\dagger}_{{{\bm{k}}}l\sigma}f_{{{\bm{k}}}l\sigma}+\sum_{{{\bm{k}}}l\sigma}\left(V^{*}_{{{\bm{k}}}l\sigma}f^{\dagger}_{{{\bm{k}}}l\sigma}c_{{{\bm{k}}}\sigma}+{\rm h.c}\right)$, where $c^{\dagger}_{{{\bm{k}}}\sigma}$ is a creation operator for $s$-electron with momentum ${{\bm{k}}}$ and spin $\sigma$ on Ce ion. $\epsilon_{{{\bm{k}}}}$ is the conduction band dispersion, which we explain in the Supplemental Material (SM) A SM . $f^{\dagger}_{{{\bm{k}}}l\Sigma}$ is a creation operator for $f$-electron with ${{\bm{k}}}$, orbital $l$ and pseudo-spin $\Sigma$. $V_{{{\bm{k}}}l\sigma}$ is the $s$-$f$ hybridization term between the nearest Ce cites. In the two-dimensional model, the pseudo-spin and $s$-electron spin are conserved ($\sigma=\Sigma$) in the $s$-$f$ mixing Tazai-HF . Using the tight-binding method Takegahara , $V_{{{\bm{k}}}l\sigma}$ is given as

[TABLE]

and $V_{{{\bm{k}}}f_{l}\downarrow}=-V_{{{\bm{k}}}f_{l}\uparrow}^{*}$. We set $A_{l}=\sqrt{{18}/{14}}$, and give a detailed explanation on $V_{{{\bm{k}}}l\sigma}$ in the SM A SM . Hereafter, we set $2|t_{ss}^{1}|=1$ as energy unit, and put $t_{sf}=0.78$, $E_{f}=-2.0$, $T=0.01$, and $\mu=-2.45$. Then, $f(s)$-electron number is $n_{f}=0.58$ ($n_{s}=0.69$).

Figure 1(a) shows the band structure of PAM. The lowest band crosses the Fermi level ($\epsilon=0$). Since $W_{D}\sim 5$eV in CeB6 ARPES1 ; ARPES2 , $2|t_{ss}^{1}|(=1)$ corresponds to $\sim 0.5$eV. The bandwidth of itinerant $f$-electron is $W_{D}^{qp}\sim|V|\sim 1$. The Fermi surfaces shown in Fig. 1(b) are composed of large ellipsoid electron pockets around X,Y points, consistently with recent ARPES studies ARPES1 ; ARPES2 .

We also introduce the Coulomb 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 normalized Coulomb interaction for Ce-ion; the maximum element of $\hat{U}^{0}$ is set to unity. The detailed explanation is given in Ref. Tazai-HF and in the SM A SM .

In the present $\Gamma_{8}$ quartet model, there are 16-type active multipole operators up to rank 3; monopole, dipole (rank 1), quadrupole (rank 2), octupole (rank 3) momenta. The table of irreducible representation (IR) for the $D_{4h}$ two-dimensional model is shown in TABLE 1 Ikeda-U . An even-rank (odd-rank) operator corresponds to an electric (magnetic) multipole operator. The $4\times 4$ matrix form of each operator, ${\hat{Q}}$, is shown in the SM B SM .

Here, we calculate the $f$-electron susceptibility. The bare irreducible susceptibility is given by $\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_{n})=({{\bm{q}}},2j\pi T)$, ${\alpha}\equiv(L,L^{\prime})$ and ${\beta}\equiv(M,M^{\prime})$. Here, ${\alpha},{\beta}$ takes $1\sim 16$, and $\hat{G}^{f}$ is the Green function without self-energy Tazai-HF . We also consider the VCs due to AL and MT terms, $\hat{X}^{\rm{AL+MT}}$, which we will explain later. Then, $f$-electron susceptibility is given as

[TABLE]

where $\hat{\phi}(q)=\hat{\chi}^{0}(q)+\hat{X}^{\rm{AL+MT}}(q)$ is irreducible susceptibility including the VCs in the $16\times 16$ matrix form.

Here, we consider the following eigen equation

[TABLE]

When the eigenvector is expressed as $\vec{w}^{\Gamma}({{\bm{q}}})=\sum_{Q\in\Gamma}Z^{Q}({{\bm{q}}})\vec{Q}$, the maximum of the eigenvalue ${\alpha}^{\Gamma}({{\bm{q}}})$ gives the Stoner factor for IR $\Gamma$, ${\alpha}^{\Gamma}=\max_{{\bm{q}}}\{{\alpha}^{\Gamma}({{\bm{q}}})\}$. Here, $\vec{Q}$ is $16\times 1$ vector defined as $(\vec{Q})_{\alpha}=(\hat{Q})_{L,L^{\prime}}$ and $Z^{Q}({{\bm{q}}})$ is a real coefficient. The $\Gamma$-channel multipole order appears when ${\alpha}^{\Gamma}\geq 1$. The inner product $(\vec{Q})^{\dagger}\vec{Q}^{\prime}$ is unity for $Q=Q^{\prime}$. It is zero when $Q$ and $Q^{\prime}$ belong to different IR, whereas it is not always zero when $Q\neq Q^{\prime}$ belong to the same IR Tazai-HF ; SM . We introduce the magnetic (electric) Stoner factor as ${\alpha}^{\rm mag(el)}=\max_{n}\{{\alpha}^{\Gamma_{n}^{-(+)}}\}$.

Using $\vec{Q}$, the multipole susceptibility is given by

[TABLE]

First, we show the numerical results by the RPA, given as $X^{\rm{AL+MT}}=0$. Figure 2 shows obtained susceptibilities $\chi^{Q}({{\bm{q}}},0)\equiv\chi^{Q,Q}({{\bm{q}}},0)$ at $u=1.08$ ($\alpha^{\rm mag}=0.9$). In the RPA, $\chi^{J_{z}}$ is the most largest. Secondly, $\chi^{T_{\nu}^{{\beta}}},\chi^{T_{\nu}^{{\alpha}}}(\nu=x,y)$ and $\chi^{T_{xyz}}$ are also enlarged. $\chi^{J_{z}}({{\bm{q}}},0)$ has peak value at ${{\bm{q}}}\approx\bm{0}$ and ${{\bm{q}}}\approx\bm{Q}\equiv(\pi,\pi)$, which is consistent with the inelastic neutron-scattering that reports strong ferromagnetic and antiferromagnetic (${{\bm{q}}}=(\pi,\pi,\pi),(\pi,\pi,0)$) fluctuations above $T_{N}$ Inosov2 ; comment2 . Therefore, the present two-dimensional PAM is reliable.

On the other hand, the RPA quadrupole susceptibility remains small. 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 quadrupole channels is much smaller than that for the dipole and octupole channels, the quadrupole susceptibilities is small within the RPA.

From now on, we introduce the VCs due to AL and MT terms. Diagrams of these VCs are shown in Fig.3 (a). For example, the AL1 term is given as

[TABLE]

where $p\equiv({{\bm{p}}},\omega_{m})$, $\bar{p}\equiv({{\bm{p}}},-\omega_{m})$, 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 is given as

[TABLE]

Other VCs are explained in the SM C SM .

Figures 3 (b) and (c) show the obtained quadrupole susceptibility by including MT- and AL-VCs. In contrast to the RPA result, the obtained $\chi^{O_{xy}}({{\bm{q}}},0)$ is strongly enhanced at ${{\bm{q}}}=\bm{Q}$ and ${{\bm{q}}}=\bm{0}$, and becomes the largest of all $\chi^{Q}$. This enhancement originates from the AL terms, whereas the MT term is very small as we show in SM C SM . The obtained $\chi^{O_{xy}}({{\bm{q}}},0)$ has the highest peak at ${{\bm{q}}}=\bm{Q}$, consistently with the antiferro-$O_{xy}$ order in CeB6. Moreover, the second highest peak of $\chi^{O_{xy}}({{\bm{q}}},0)$ at ${{\bm{q}}}=\bm{0}$ explains the softening of shear modulus $C_{44}$ in CeB6 Goto . We show other quadrupole susceptibilities in the SM C SM . To summarize, the obtained strong enhancements of $\chi^{O_{xy}}({{\bm{q}}},0)$ and $\chi^{J_{z}}({{\bm{q}}},0)$ at both ${{\bm{q}}}=\bm{Q}$ and ${{\bm{q}}}=\bm{0}$ reproduce the key experimental results of CeB6.

Next, we explain that the $O_{xy}$ quadrupole order is derived from the interference between magnetic multipole fluctuations. For this purpose, we analyze the total AL term $\hat{X}\equiv\hat{X}^{\rm AL1}+\hat{X}^{\rm AL2}$ for $O_{xy}$-channel defined as

[TABLE]

where $\hat{X}\equiv\hat{X}^{\rm AL1}+\hat{X}^{\rm AL2}$. The Stoner factor for $O_{xy}(=\Gamma_{4}^{+})$ channel is proportional to $uU^{O_{xy}}_{0}\phi_{O_{xy}}(q)$, where $\phi_{O_{xy}}(q)\equiv(\vec{O}_{xy})^{\dagger}\hat{\phi}(q)\vec{O}_{xy}$. Therefore, $X_{O_{xy}}(q)\ (>0)$ works as enhancement factor of $O_{xy}$ susceptibility.

By following Ref. Tazai-HF , we expand $\hat{V}(q)$ on the basis of multipole operator as

[TABLE]

where the real coefficient $v^{QQ^{\prime}}_{q}$ is uniquely determined Tazai-HF . From Eq.(7), (9) and (10), the AL1 term due to $(Q,Q^{\prime})$-channel fluctuations is given as

[TABLE]

where $v^{Q}\equiv v^{QQ}$ and $\Lambda_{q,p}^{O_{xy}QQ^{\prime}}$ is defined as

[TABLE]

The diagrammatic expression of Eq. (11) is shown in Fig.4(a). Figure 4(b) shows the ${{\bm{q}}}$-dependence of $X_{O_{xy}}^{QQ^{\prime}}({{\bm{q}}},0)$ at $u=0.91$. We find that the $(Q,Q^{\prime})=(T_{x}^{{\alpha}},T_{y}^{{\alpha}})$, $(J_{z},T_{xyz})$, $(T_{x}^{{\beta}},T_{y}^{{\beta}})$ channels give the dominant contributions. Other terms not shown in Fig.4(b) give negligible contribution.

Figure 4(c) presents the quantum process of $O_{xy}$ quadrupole order driven by the interference between $(T_{x},T_{y})$ fluctuations, which corresponds to $\Lambda^{O_{zx}T_{x}T_{y}}$ in Fig.4(a). This process is realized when $\Lambda^{O_{zx}QQ^{\prime}}\sim{\rm Tr}\{\hat{O}_{xy}\cdot\hat{Q}\cdot\hat{Q}^{\prime}\}\neq 0$. Since $\Lambda^{QTT^{\prime}}=0$ for odd-rank $Q$, the AL-VC is unimportant for $\chi^{J}$ and $\chi^{T}$ Yamakawa-FeSe .

Next, the ${{\bm{q}}}$-dependence of the AL-VC is given as $X_{O_{xy}}^{T_{x}T_{y}}({{\bm{q}}})\propto\sum_{{{\bm{p}}}}\chi^{T_{x}}({{\bm{p}}})\chi^{T_{y}}({{\bm{q}}}-{{\bm{p}}})$, which becomes large at ${{\bm{q}}}=\bm{Q}$ and ${{\bm{q}}}=\bm{0}$ since $\chi^{T_{\mu}}({{\bm{p}}})$ has large peaks at ${{\bm{p}}}\sim\bm{Q},\bm{0}$ shown in Fig. 2. Thus, antiferro-quadrupole order in CeB6 originates from the interference between ferro- and antiferro-magnetic multipole fluctuations.

Finally, we discuss the field-induced octupole order, which has been studied intensively as a main issue of CeB6 Shiina1 ; Shiina2 ; Shiina3 ; Shiina4 . The Zeeman term under the magnetic field along $z$-axis is given as $\hat{H}_{Z}=h_{z}\sum_{L,M}(\hat{J}_{z})_{L,M}f^{\dagger}_{{{\bm{k}}}L}f_{{{\bm{k}}}M}$. When $h_{z}\neq 0$, both $O_{xy}$ and $T_{xyz}$ belong to the same IR $\Gamma_{4}$ shown in TABLE 1 Shiina1 . Therefore, large quadrupole-octupole susceptibility $\chi^{O_{xy},T_{xyz}}({{\bm{q}}},0)$ is induced in proportion to $h_{z}$. To verify this, we solve the eigen equation (4) for the IR $\Gamma_{4}$ under $h_{z}$, at the fixed magnetic Stoner factor in the RPA ${\alpha}^{\rm mag}=0.8$ comment ; Sakurazawa .

Figures 5(a) and (b) show the obtained eigenvector $\vec{w}^{\Gamma_{4}}({{\bm{q}}})=Z^{O_{xy}}({{\bm{q}}})\vec{O}_{xy}+Z^{T_{xyz}}({{\bm{q}}})\vec{T}_{xyz}$ ($|\vec{w}^{\Gamma_{4}}|^{2}=1$) and the Stoner factor ${\alpha}^{\Gamma_{4}}$ at ${{\bm{q}}}=\bm{Q}$, respectively, as functions of $h_{z}$. Here, ${\alpha}^{\Gamma_{4}}$ is the largest Stoner factor. The increment of ${\alpha}^{\Gamma_{4}}$ under $h_{z}$ is consistent with the field-enhancement of $T_{Q}$ in CeB6. (In contrast, $T_{N}$ will be suppressed by large $O_{xz}$ moment.) Also, $Z^{T_{xyz}}$ increases linearly in $h_{z}$, due to the interference process under $h_{z}$ shown in the inset of Fig. 5(b). $Z^{T_{xyz}}$ becomes comparable to $Z^{O_{xy}}$ under small magnetic field $h_{z}\lesssim 0.03\ll W_{D}^{qp}/10$. Since the ratio of the ordered momenta at $T_{Q}$ is $M^{T_{xyz}}/M^{O_{xy}}=Z^{T_{xyz}}/Z^{O_{xy}}$, field-induced antiferro-$T_{xyz}$ order is naturally explained.

In summary, we developed multipole fluctuation theory by focusing on the AL-type VCs in HF systems, and applied the theory to the multipole order physics in CeB6. Both ferro- and antiferro-magnetic multipole fluctuations emerge in CeB6 due to the nesting of Fermi surfaces, consistently with neutron experiments. Then, antiferro-$O_{xy}$ order in CeB6 at $T_{Q}\ (>T_{N})$ is derived from the interference between different magnetic multipole fluctuations, which is depicted in Fig. 4 (c). We also explained the field-induced octupole order, which is a central issue of CeB6. The discovered inter-multipole coupling mechanism will be significant in various HF systems, such as quadrupole ordering system Pr$T_{2}$Zn20 ($T$ = Rh and Ir) Oni-Pr and Pr$T_{2}$Al20 ($T$=V,Ti) Naka-Pr . Although the analysis of AL-VC in three-dimensional PAM is very difficult, it is an important future problem.

We stress that the on-site quadrupole ($O_{xy}$) interaction on Ce-ion is about 60% of dipole ($J_{\mu}$) one as shown in TABLE 2. Therefore, quadrupole order cannot appear within the mean-field theory. In contrast, in the localized RKKY model, quadrupole interaction is as large as the dipole interaction Shiina1 ; Shiina4 . Such discrepancy between itinerant picture and localized one, which is an important problem in HF systems, is partially resolved by considering the VCs as we discussed here.

I.1 A: model Hamiltonian

Here, we present detailed explanation for the model Hamiltonian. In CeB6, the conduction band is composed of $5d$ electrons on Ce-ions, Here, to simplify the model Hamiltonian, we introduce the conduction band made of $s$ electrons. The realistic tight-binding model of conduction band of CeB6 is given in Ref. S-ARPES2 . In the present study, we slightly modify the model in Ref. S-ARPES2 and put $k_{z}=0$, in order to reproduce the experimental Fermi surfaces of CeB6 on the $k_{x}$-$k_{y}$ plane after $s$-$f$ hybridization. The present two-dimensional tight-binding model for conduction band is given as

[TABLE]

where $t_{ss}^{i}$ is the $i$-th nearest $s$-$s$ hopping integral. We set $(t_{ss}^{1},t_{ss}^{2},t_{ss}^{3},t_{ss}^{4},t_{ss}^{5})=(-0.5,-0.889,0.292,-0.229,0.687)$, and $E_{0}=1.33$.

Next, we explain the hybridization term. Based on the Slater-Koster tight-binding method, the $s$-$f$ hybridization between the nearest Ce-sites is

[TABLE]

and $V_{{{\bm{k}}}f_{l}\downarrow}=-V_{{{\bm{k}}}f_{l}\uparrow}^{*}$. Here, $t_{sf}=(sf\sigma)$, and $A_{1}=\sqrt{{18}/{14}}$ and $A_{2}=\sqrt{{3}/{7}}$. Since $A_{1}>A_{2}$, the relation $D_{f_{1}}(0)>D_{f_{2}}(0)$ holds in the present two-dimensional PAM, where $D_{f_{l}}(0)$ is the $f_{l}$-electron density-of-states at Fermi level. However, $D_{f_{1}}(0)=D_{f_{2}}(0)$ holds in the cubic model, since the $s$-$f$ hybridization along $z$-axis is larger for $f_{2}$-electron. To escape from the artifact of two-dimensionality, we put $A_{1}=A_{2}=\sqrt{{18}/{14}}$ in the present study.

In the present $\Gamma_{8}$ model, the relation $V_{{{\bm{k}}}f_{1}{\sigma}}\propto V_{{{\bm{k}}}f_{2}{\sigma}}^{*}$ holds as shown in Eq. (S2). In contrast, in the $\Gamma_{7}^{(1)}$-$\Gamma_{7}^{(2)}$ model for CeCu2Si2 used in Ref. S-Tazai-HF , the relation $V_{{{\bm{k}}}f_{1}{\sigma}}\propto V_{{{\bm{k}}}f_{2}{\sigma}}$ holds.

Finally, we explain the Coulomb interaction in $f$-electrons, which is derived from Slater-Condon parameter $F^{p}$ S-Tazai-HF . We set $(F^{0},F^{2},F^{4},F^{6})=(5.3,9.09,6.927,4.756)$ in unit eV by referring Ref.S-F0F2F4F6 . The derived Coulomb interaction is about 6eV. If we use the such large Coulomb interaction in the RPA, the magnetic order appears since the self-energy is dropped in the RPA. Therefore, we introduce the following Coulomb interaction term:

[TABLE]

where $L=(l,\sigma)$ and $M=(m,\rho)$. $u$ is the interaction model parameter, and $\hat{U}^{0}$ is the normalized Coulomb interaction introduced in Ref. S-Tazai-HF . That is, the maximum element of $\hat{U}^{0}$ is normalized to unity.

I.2 B: multipole-operator

Here, we list the pseudo-spin representation of the multipole operators in TABLE 1, which was first introduced in Ref.S-Shiina1 . An even-rank (odd-rank) operator corresponds to an electric (magnetic) multipole operator. Each multipole operator of rank $k$ are composed of $4\times 4$ tensor $J^{(k)}_{q}(q=-k\sim k)$ S-Shiina1 ; S-Springer which is given by $[J_{\pm},J^{(k)}_{q}]=\sqrt{(k\mp q)(k\pm q+1)}J^{(k)}_{q\pm 1}$ $J_{k}^{(k)}=(-1)^{k}\sqrt{(2k-1)!!/(2k)!!}J_{+}^{k}$. The multipole operators $\hat{Q}$ is given by the linear combination of $J^{(k)}_{q}$. The $4\times 4$ matrix form of each electric (odd-rank) multipole operators is given by S-Shiina1

[TABLE]

The $4\times 4$ matrix form of each magnetic (odd-rank) multipole operators is given by S-Shiina1

[TABLE]

In the main text, we use the normalized multipole matrix introduced as follows:

[TABLE]

Then, the normalized $\hat{Q}$ satisfies the condition $\sum_{L,M}|Q_{L,M}|^{2}=1$.

I.3 C: multipole fluctuations

In the main text, we explain the analytic expression only for AL1 term. The expression for the AL2 term is given as

[TABLE]

where

[TABLE]

The expression for the MT term is

[TABLE]

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

In the main text, we perform the numerical study of multipole susceptibilities by considering both MT- and AL-VCs, and showed that $O_{xy}$ octupole susceptibility is strongly enlarged by the AL-VCs. Here, we show all the quadrupole susceptibilities obtained by the present study in Fig. S1. In the cubic model, $\chi^{Q}({{\bm{q}}},0)$ with $Q=O_{xy},O_{zx},O_{yz}$ should equally develop. In the present two-dimensional model, however, only $O_{xy}$-fluctuation strongly develops. The reason is that ($T_{x},T_{y}$) fluctuations are much larger than $T_{z}$ fluctuations in the RPA, due to the violation of cubic symmetry. Since $O_{\mu\nu}$ quadrupole susceptibility is magnified by ($T_{\mu},T_{\nu})$ fluctuations ($\mu,\nu=x,y,z$) due to the AL-VC, $\chi^{O_{xy}}({{\bm{q}}},0)$ is the largest in the present model.

As we show in TABLE 2, the Coulomb interaction $U^{Q}_{0}$ for $Q=O_{xy/yz/zx}$ is much larger than that for $Q=O_{20/22}$. For this reason, it is difficult to expect that $Q=O_{20/22}$ quadrupole susceptibility becomes larger than $Q=O_{xy}$ one, even if the AL-VCs are considered. Thus, the relation $\chi^{O_{xy}}({{\bm{q}}},0)>\chi^{O_{20/22}}({{\bm{q}}},0)$ should hold even in cubic systems.

Next, we calculate the susceptibility with AL-VC (MT-VC), $\chi^{O_{xy}}_{\rm{AL(MT)}}({{\bm{q}}},0)$, given by $\hat{\phi}(q)=\hat{\chi}^{0}(q)+\hat{X}^{\rm{AL(MT)}}(q)$. Figure S2 shows the obtained $\chi^{O_{xy}}_{\rm{AL}}({{\bm{q}}},0)$ and $\chi^{O_{xy}}_{\rm{MT}}({{\bm{q}}},0)$ as functions of $u$. $\chi^{O_{xy}}_{\rm{AL}}({{\bm{q}}},0)$ strongly increases with $u$, similarly to $\chi^{O_{xy}}({{\bm{q}}},0)$ with AL+MT terms shown in Fig. 3 (c) in the main text. In contrast, $\chi^{O_{xy}}_{\rm{MT}}({{\bm{q}}},0)$ remains small and comparable to the RPA result in Fig. 3 (c). Therefore, it is verified that the enhancement of $O_{xy}$ quadrupole fluctuations originates from the AL-VC, whereas the MT-VC is very small.

To understand this result analytically, we analyze the AL and MT terms for the electric multipole channel given by the following magnetic multipole susceptibility

[TABLE]

where $\xi^{2}\propto(T-T_{0})^{-1}$ and ${\omega}_{\rm mag}\propto\xi^{-2}$. $\xi$ is the correlation length. Then, in two-dimensional systems at a fixed $T$, AL-VC and MT-VC given in Eqs. (S8)-(S9) are scaled as $X^{\rm AL}(\bm{0},0)\sim\sum_{{\bm{p}}}\{\chi^{\rm mag}({{\bm{p}}},0)\}^{2}\sim\xi^{2}$ and $X^{\rm MT}(\bm{0},0)\sim\sum_{{\bm{p}}}\chi^{\rm mag}({{\bm{p}}},0)\sim\log\xi$, respectively. Therefore, the AL term dominates over the MT term when $\xi\gg 1$ S-Onari . The significance of the AL terms near the magnetic criticality is verified by the functional-renormalization-group (fRG) study S-RG1 ; S-RG2 ; S-RG3 .

In $d$-dimensional system, the AL term is proportional to $\max\{\xi^{4-d},1\}$. This fact means that the non-locality of irreducible AL diagram is significant near the magneitc criticality $(\xi\gg 1)$.