Derivation of RKKY Interaction between Multipole Moments in CeB$_6$ by the Effective Wannier Model based on the Bandstructure Calculation

TL;DR

This study calculates the RKKY interactions among multipole moments in CeB$_6$ using a detailed Wannier model derived from bandstructure calculations, explaining various ordered phases and providing a method applicable to similar materials.

Contribution

The paper introduces a realistic Wannier-based approach to compute RKKY interactions for all multipole moments in CeB$_6$, linking bandstructure data to magnetic ordering phenomena.

Findings

Enhanced quadrupole and octupole couplings at specific wavevectors explain phase II ordering.

Large couplings for certain octupoles suggest possible phase IV ordering.

Intersite RKKY couplings show diverse long-range behaviors depending on multipole type.

Abstract

We have investigated the electronic states of CeB$_6$ and have directly calculated the RKKY interaction on the basis of the 74-orbital effective Wannier model which includes 14 Ce-$f$ orbitals and 60 conduction ($c$) orbitals of Ce-$d,s$ and B-$p,s$ derived from the density-functional theory bandstructure calculation. By using not only the $c$-band dispersion but also the $f$-$c$ mixing matrix elements of the Wannier model, the realistic couplings for all 15 active multipole moments in $\Gamma_8$ quartet subspace are obtained in the wavevector $q$-space and real-space. Both of the $\Gamma_{5g}$ quadrupoles $(O_{yz},O_{zx},O_{xy})$ and the $\Gamma_{2u}$ octupole $T_{xyz}$ couplings are maximally enhanced with $q=(\pi,\pi,\pi)$ which naturally explains the phase II of the antiferro-quadrupolar ordering at $T_{Q}=3.2$ K, and are also enhanced with $q=(0,0,0)$ corresponding to the elastic softening of $C_{44}$. Also the couplings of the $\Gamma_{5u}$ octupoles $T_{x}^{\beta}$, $T_{y}^{\beta}$ and $T_{z}^{\beta}$ are quite large for $q=(\pi,0,0)$, $(0,\pi,0)$ and $(0,0,\pi)$, which yields the antiferro-octupolar ordering of a possible candidate for phase IV of Ce$_{x}$La$_{1-x}$B$_6$. The intersite vector dependence of the RKKY couplings exhibit different long-range, oscillating, isotropic and anisotropic behaviors depending on the types of the multipole moments. The present approach enables us to provide the information about the possible multipole ordering in an unbiased way and is easily available for other localized $f$ electron materials once the $c$ states and $f$-$c$ mixing elements are given from the bandstructure calculation.

Figures

Tables

Table 1. Table 1: The irreducible representations (IRRs) and notations for the active multipole moments in Γ 8 subscript Γ 8 \Gamma_{8} subspace [ 3 ] where 𝑱 ​ ( 𝑻 α , β ) 𝑱 superscript 𝑻 𝛼 𝛽 \bm{J}~{}(\bm{T}^{\alpha,\beta}) is the dipole (octupole), 𝑱 = ( J x , J y , J z ) 𝑱 subscript 𝐽 𝑥 subscript 𝐽 𝑦 subscript 𝐽 𝑧 \bm{J}=(J_{x},J_{y},J_{z}) , 𝑻 α ​ ( β ) = ( T x α ​ ( β ) , T y α ​ ( β ) , T z α ​ ( β ) ) superscript 𝑻 𝛼 𝛽 superscript subscript 𝑇 𝑥 𝛼 𝛽 superscript subscript 𝑇 𝑦 𝛼 𝛽 superscript subscript 𝑇 𝑧 𝛼 𝛽 \bm{T}^{\alpha(\beta)}=(T_{x}^{\alpha(\beta)},T_{y}^{\alpha(\beta)},T_{z}^{\alpha(\beta)}) , η ± = − ( τ z ∓ 3 ​ τ x ) / 2 superscript 𝜂 plus-or-minus minus-or-plus superscript 𝜏 𝑧 3 superscript 𝜏 𝑥 2 \eta^{\pm}=-(\tau^{z}\mp\sqrt{3}\tau^{x})/2 , ζ ± = − 1 2 ​ ( τ x ± 3 ​ τ z ) superscript 𝜁 plus-or-minus 1 2 plus-or-minus superscript 𝜏 𝑥 3 superscript 𝜏 𝑧 \zeta^{\pm}=-\frac{1}{2}(\tau^{x}\pm\sqrt{3}\tau^{z}) , and g ​ ( u ) 𝑔 𝑢 g~{}(u) means even (odd) time-reversal symmetry. In this paper we call 2 9 ​ 5 ​ T x ​ y ​ z 2 9 5 subscript 𝑇 𝑥 𝑦 𝑧 \frac{2}{9\sqrt{5}}T_{xyz} , 1 4 ​ ( O u , O v ) 1 4 subscript 𝑂 𝑢 subscript 𝑂 𝑣 \frac{1}{4}(O_{u},O_{v}) and 1 3 ​ 5 ​ 𝑻 β 1 3 5 superscript 𝑻 𝛽 \frac{1}{3\sqrt{5}}\bm{T}^{\beta} just as T x ​ y ​ z subscript 𝑇 𝑥 𝑦 𝑧 T_{xyz} , ( O u , O v ) subscript 𝑂 𝑢 subscript 𝑂 𝑣 (O_{u},O_{v}) and 𝑻 β superscript 𝑻 𝛽 \bm{T}^{\beta} , but all the multipole operators are normalized.

IRR [dimension]	vector	pseudospin	multipole
${\Gamma }_{2u}$ [1]	$\xi$	${\tau }^y$	$\frac{2}{9\sqrt{5}}{T}_{xyz}$
${\Gamma }_{3g}$ [2]	${𝝉}^{\prime }$	$({\tau }^z,{\tau }^x)$	$\frac{1}{4}(O_u,O_v)$
${\Gamma }_{5g}$ [3]	$𝝁$	$({\tau }^y{\sigma }^x,{\tau }^y{\sigma }^y,{\tau }^y{\sigma }^z)$	$(O_{yz},O_{zx},O_{xy})$
${\Gamma }_{4u}^{\left(1\right)}$ [3]	$𝝈$	$({\sigma }^x,{\sigma }^y,{\sigma }^z)$	$\frac{14}{15}𝑱-\frac{4}{45}{𝑻}^{\alpha }$
${\Gamma }_{4u}^{\left(2\right)}$ [3]	$𝜼$	$({\eta }^{+}{\sigma }^x,{\eta }^{-}{\sigma }^y,{\tau }^z{\sigma }^z)$	$-\frac{2}{15}𝑱+\frac{7}{45}{𝑻}^{\alpha }$
${\Gamma }_{5u}$ [3]	$𝜻$	$({\zeta }^{+}{\sigma }^x,{\zeta }^{-}{\sigma }^y,{\tau }^x{\sigma }^z)$	$\frac{1}{3\sqrt{5}}{𝑻}^{\beta }$

Table 2. Table 2: The obtained possible multipole modes upto the 10th largest coupling together with the corresponding moment types, wavevectors, maximum values and ratios to the largest value of Γ 5 ​ g subscript Γ 5 𝑔 \Gamma_{5g} - ( π , π , π ) 𝜋 𝜋 𝜋 (\pi,\pi,\pi) , where the Γ 3 ​ g subscript Γ 3 𝑔 \Gamma_{3g} moments O y 2 − z 2 subscript 𝑂 superscript 𝑦 2 superscript 𝑧 2 O_{y^{2}-z^{2}} , O z 2 − x 2 subscript 𝑂 superscript 𝑧 2 superscript 𝑥 2 O_{z^{2}-x^{2}} , O 3 ​ x 2 − r 2 subscript 𝑂 3 superscript 𝑥 2 superscript 𝑟 2 O_{3x^{2}-r^{2}} and O 3 ​ y 2 − r 2 subscript 𝑂 3 superscript 𝑦 2 superscript 𝑟 2 O_{3y^{2}-r^{2}} are described by the linear combinations of O u ( = O 3 ​ z 2 − r 2 ) annotated subscript 𝑂 𝑢 absent subscript 𝑂 3 superscript 𝑧 2 superscript 𝑟 2 O_{u}(=O_{3z^{2}-r^{2}}) and O v ( = O x 2 − y 2 ) annotated subscript 𝑂 𝑣 absent subscript 𝑂 superscript 𝑥 2 superscript 𝑦 2 O_{v}(=O_{x^{2}-y^{2}}) .

rank	IRR	multipole	wavevector	value [meV]	ratio
1	${\Gamma }_{5g}$	$(O_{yz},O_{zx},O_{xy})$	$(\pi ,\pi ,\pi )$	17.26	1.00
\hdashline 2	${\Gamma }_{2u}$	$T_{xyz}$	$(\pi ,\pi ,\pi )$	16.56	0.96
3	${\Gamma }_{5u}$	${\zeta }^x$ ${\zeta }^y$ ${\zeta }^z$	$(\pi ,0,0)$ $(0,\pi ,0)$ $(0,0,\pi )$	14.48	0.84
\hdashline 4	${\Gamma }_{3g}$	$O_{y^2-z^2}$ $O_{z^2-x^2}$ $O_{x^2-y^2}$	$(\pi ,0,0)$ $(0,\pi ,0)$ $(0,0,\pi )$	14.08	0.82
5	${\Gamma }_{2u}$	$T_{xyz}$	$(0,0,0)$	12.43	0.72
\hdashline 6	${\Gamma }_{5g}$	$(O_{yz},O_{zx},O_{xy})$	$(0,0,0)$	11.69	0.68
7	${\Gamma }_{4u}^{\left(1\right)}$	${\sigma }^x$ ${\sigma }^y$ ${\sigma }^z$	$(0,\pi ,0),(0,0,\pi )$ $(0,0,\pi ),(\pi ,0,0)$ $(\pi ,0,0),(0,\pi ,0)$	8.41	0.49
8	${\Gamma }_{4u}^{\left(1\right)}$	${\sigma }^x$ ${\sigma }^y$ ${\sigma }^z$	$(\pi ,0,0)$ $(0,\pi ,0)$ $(0,0,\pi )$	8.11	0.47
9	${\Gamma }_{4u}^{\left(2\right)}$	${\eta }^x$ ${\eta }^y$ ${\eta }^z$	$(0,\pi ,0),(0,0,\pi )$ $(0,0,\pi ),(\pi ,0,0)$ $(\pi ,0,0),(0,\pi ,0)$	7.87	0.46
\hdashline 10	${\Gamma }_{3g}$	$O_{3x^2-r^2}$ $O_{3y^2-r^2}$ $O_{3z^2-r2}$	$(0,\pi ,0),(0,0,\pi )$ $(0,0,\pi ),(\pi ,0,0)$ $(\pi ,0,0),(0,\pi ,0)$	7.78	0.45

Equations

$$H_{\mathrm{TB}}$$

$$+\sum _{ij}\sum _{m\ell }\left(V_{im,j\ell }{f}_{im}^{†}{c}_{j\ell }+h.c.\right),$$

$$H_{\mathrm{TB}}$$

$$+\sum _{\mathbf{k}}\sum _{m\ell }\left(V_{\mathbf{k}m\ell }{f}_{\mathbf{k}m}^{†}{c}_{\mathbf{k}\ell }+h.c.\right)=\sum _{\mathbf{k}s}{\epsilon }_{\mathbf{k}s}{a}_{\mathbf{k}s}^{†}{a}_{\mathbf{k}s},$$

$${\Sigma }_m^f(\epsilon )=\Delta {\epsilon }_m^f+\left(1-\frac{1}{Z_m^f}\right)\epsilon -i{\gamma }_m{\epsilon }^2+O({\epsilon }^3),$$

$$\displaystyle\Delta\varepsilon_{m}^{f}={\rm Re}\Sigma_{m}^{f}(0),~{}~{}Z_{m}^{f}=\left(1-\frac{d}{d\varepsilon}{\rm Re}\Sigma_{m}^{f}(\varepsilon)\Bigr{|}_{\varepsilon=0}\right)^{-1},$$

$$H_{\mathrm{RTB}}$$

$$+\sum _{ij}\sum _{m\ell }\left({\tilde{V}}_{im,j\ell }{f}_{im}^{†}{c}_{j\ell }+h.c.\right)$$

$$\displaystyle\tilde{h}_{im,jm^{\prime}{}}^{ff}=\left\{\begin{array}[]{cc}\varepsilon_{m}^{f}+\Delta\varepsilon_{m}^{f}&i=j,~{}m=m^{\prime}{}\\ \sqrt{Z_{m}^{f}Z_{m^{\prime}{}}^{f}}h_{im,jm^{\prime}{}}^{ff}&i\neq j\\ \end{array}\right.$$

$${\tilde{V}}_{im,j\ell }=\sqrt{Z_m^f}{V}_{im,j\ell }$$

$$m_{\mathrm{spin}}=-{\mu }_{\mathrm{B}}\frac{1}{N}\sum _{\mathbf{k}s}\sum _{m{m}^{\prime }}{\left({\sigma }^z\right)}_{m{m}^{\prime }}{u}_{\mathbf{k}sm}{u}_{\mathbf{k}s{m}^{\prime }}^{\ast }f\left({\epsilon }_{\mathbf{k}s}\right),$$

$$m_{\mathrm{orb}}=-{\mu }_{\mathrm{B}}\frac{1}{N}\sum _{\mathbf{k}s}\sum _{m{m}^{\prime }}{\left({\ell }^z\right)}_{m{m}^{\prime }}{u}_{\mathbf{k}sm}{u}_{\mathbf{k}s{m}^{\prime }}^{\ast }f\left({\epsilon }_{\mathbf{k}s}\right),$$

$${\chi }_{O_{{}_{\Gamma }}}(\mathbf{q})=\sum _{m_1{m}_2}\sum _{m_3{m}_4}{O}_{m_1{m}_2}^{\Gamma }{O}_{m_4{m}_3}^{\Gamma }{\chi }_{m_1{m}_2{m}_3{m}_4}(\mathbf{q}),$$

$${\chi }_{m_1{m}_2{m}_3{m}_4}(\mathbf{q})=\frac{1}{N}\sum _{\mathbf{k}s{s}^{\prime }}{u}_{\mathbf{k}s{m}_3}^{\ast }{u}_{\mathbf{k}s{m}_1}{u}_{\mathbf{k}+\mathbf{q}{s}^{\prime }{m}_2}^{\ast }{u}_{\mathbf{k}+\mathbf{q}{s}^{\prime }{m}_4}$$

$$\times \frac{f({\epsilon }_{\mathbf{k}+\mathbf{q}{s}^{\prime }})-f({\epsilon }_{\mathbf{k}s})}{{\epsilon }_{\mathbf{k}s}-{\epsilon }_{\mathbf{k}+\mathbf{q}{s}^{\prime }}},$$

$$H_{\mathrm{MKL}}$$

$$+\sum _i\sum _{m{m}^{\prime }}\sum _{\mathbf{k}{\mathbf{k}}^{\prime }}\sum _{\ell {\ell }^{\prime }}{J}_{im{m}^{\prime }}^{\mathbf{k}\ell ,{\mathbf{k}}^{\prime }{\ell }^{\prime }}{f}_{im}^{†}{f}_{i{m}^{\prime }}{c}_{\mathbf{k}\ell }^{†}{c}_{{\mathbf{k}}^{\prime }{\ell }^{\prime }},$$

$$J_{im{m}^{\prime }}^{\mathbf{k}\ell ,{\mathbf{k}}^{\prime }{\ell }^{\prime }}=\frac{2}{N}\frac{V_{\mathbf{k}m\ell }{V}_{{\mathbf{k}}^{\prime }{m}^{\prime }{\ell }^{\prime }}^{\ast }}{\mu -{\epsilon }_{{\Gamma }_8}^f}{e}^{-i(\mathbf{k}-{\mathbf{k}}^{\prime })\cdot {\mathbf{R}}_i}$$

$$H_{\mathrm{RKKY}}=-\sum _{⟨ij⟩}\sum _{m_1{m}_2}\sum _{m_3{m}_4}{K}_{m_1{m}_2{m}_3{m}_4}({\mathbf{R}}_{ij})f_{i{m}_1}^{†}{f}_{i{m}_2}{f}_{j{m}_4}^{†}{f}_{j{m}_3},$$

$$K_{m_1{m}_2{m}_3{m}_4}({\mathbf{R}}_{ij})=\frac{1}{N}\sum _{\mathbf{q}}{K}_{m_1{m}_2{m}_3{m}_4}(\mathbf{q})e^{i\mathbf{q}\cdot ({\mathbf{R}}_i-{\mathbf{R}}_j)},$$

$$K_{m_1{m}_2{m}_3{m}_4}(\mathbf{q})=\frac{1}{N}\sum _{\mathbf{k}s{s}^{\prime }}\sum _{{\ell }_1{\ell }_2}\sum _{{\ell }_3{\ell }_4}\frac{V_{\mathbf{k}{m}_3{\ell }_3}^{\ast }{V}_{\mathbf{k}{m}_1{\ell }_1}{V}_{\mathbf{k}+\mathbf{q}{m}_2{\ell }_2}^{\ast }{V}_{\mathbf{k}+\mathbf{q}{m}_4{\ell }_4}}{(\mu -{\epsilon }_{{\Gamma }_8}^f{)}^2}$$

$$\times u_{\mathbf{k}s{\ell }_3}^{c\ast }{u}_{\mathbf{k}s{\ell }_1}^c{u}_{\mathbf{k}+\mathbf{q}{s}^{\prime }{\ell }_2}^{c\ast }{u}_{\mathbf{k}+\mathbf{q}{s}^{\prime }{\ell }_4}^c\frac{f({\epsilon }_{\mathbf{k}+\mathbf{q}{s}^{\prime }})-f({\epsilon }_{\mathbf{k}s})}{{\epsilon }_{\mathbf{k}s}-{\epsilon }_{\mathbf{k}+\mathbf{q}{s}^{\prime }}},$$

$$v_{\mathbf{k}s}^{m{m}^{\prime }}=\sum _{\ell {\ell }^{\prime }}{V}_{\mathbf{k}m\ell }^{\ast }{V}_{\mathbf{k}{m}^{\prime }{\ell }^{\prime }}{u}_{\mathbf{k}s\ell }^{c\ast }{u}_{\mathbf{k}s{\ell }^{\prime }}^c,$$

$$K_{m_1{m}_2{m}_3{m}_4}(\mathbf{q})=\frac{1}{N}\sum _{\mathbf{k}s{s}^{\prime }}\frac{v_{\mathbf{k}s}^{m_3{m}_1}{v}_{\mathbf{k}+\mathbf{q}{s}^{\prime }}^{m_2{m}_4}}{(\mu -{\epsilon }_{{\Gamma }_8}^f{)}^2}\frac{f({\epsilon }_{\mathbf{k}+\mathbf{q}{s}^{\prime }})-f({\epsilon }_{\mathbf{k}s})}{{\epsilon }_{\mathbf{k}s}-{\epsilon }_{\mathbf{k}+\mathbf{q}{s}^{\prime }}}.$$

$$H_{\mathrm{RKKY}}^{\mathrm{MF}}=-\sum _{\mathbf{q}}\sum _{\Gamma }{\overline{K}}_{O_{\Gamma }}(\mathbf{q}){\overline{O}}_{\Gamma }(\mathbf{q}){\overline{O}}_{\Gamma }(-\mathbf{q}),$$

$$K_{O_{\Gamma }}(\mathbf{q})=\sum _{m_1{m}_2}\sum _{m_3{m}_4}{O}_{m_1{m}_2}^{\Gamma }{O}_{m_4{m}_3}^{\Gamma }{K}_{m_1{m}_2{m}_3{m}_4}(\mathbf{q}),$$

$${\overline{O}}_{\Gamma }(\mathbf{q})=\frac{1}{N}\sum _i{\overline{O}}_{\Gamma }({\mathbf{R}}_i)e^{i\mathbf{q}\cdot {\mathbf{R}}_i},$$

$${\chi }_{O_{\Gamma }}^{\mathrm{MF}}(\mathbf{q})=\frac{{\chi }_{O_{\Gamma }}(\mathbf{q})}{1-{\chi }_{O_{\Gamma }}(\mathbf{q}){\overline{K}}_{O_{\Gamma }}(\mathbf{q})},$$

$$K_{O_{\Gamma }}({\mathbf{R}}_{ij})=\sum _{m_1{m}_2}\sum _{m_3{m}_4}{O}_{m_1{m}_2}^{\Gamma }{O}_{m_4{m}_3}^{\Gamma }{K}_{m_1{m}_2{m}_3{m}_4}({\mathbf{R}}_{ij}),$$

Full text

1 Introduction

CeB6 has been known as a typical and remarkable compound exhibiting a rich phase diagram of the multipole orderings [1, 2, 3, 4] and extensively studied experimentally [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31] and theoretically[32, 33, 34, 35, 36, 37, 38, 39, 40, 41]. Due to the large spin-orbit coupling (SOC) and the cubic crystalline electric field (CEF), the ground state of 4$f^{1}$ in Ce3+ ion is the $\Gamma_{8}$ quartet separated from the excited $\Gamma_{7}$ doublet by 540 K[8], and has a inherently the degrees of freedom of 15 active multipole moments as shown in Table 1.

Up to now, three phases exist in temperature $T$ and external magnetic field $H$ plane of CeB6. Normal phase (phase I) from a room temperature down to a few K with $H=0$ is a typical Kondo lattice metal with a highly-enhanced specific-heat coefficient[9, 10] $C/T=250$ mJ/mol$\cdot$K2. With decreasing $T$, phase II emerges at a critical temperature $T_{Q}=3.2$ K with the ordering wavevector $\bm{q}=(\pi,\pi,\pi)$ and is confirmed by the antiferro-quadrupolar (AFQ) ordering of the $\Gamma_{5g}$ quadrupoles $(O_{yz},O_{zx},O_{xy})$. The ordering tendency of the $\Gamma_{5g}$ quadrupole moment is supported from the elastic-softening of $C_{44}$ at low temperature[15, 16, 17, 18]. Interestingly, $T_{Q}$ increases with increasing the applied field $H$, where the $\Gamma_{2u}$ octupole $T_{xyz}$ moment is induced by $H$ in addition to the $\Gamma_{5g}$ quadrupoles, which is well understood by the analysis of NMR[21]. The phase III is a antiferro-magnetic (AFM) ordering of $\Gamma_{4u}$ magnetic moments $(\sigma^{x},\sigma^{y},\sigma^{z})$ at $T_{N}=2.3$ K with the double-$\bm{q}$-structure of $\bm{Q}_{1}=(\tfrac{\pi}{2},\tfrac{\pi}{2},0)$ and $\bm{Q}_{2}=(\tfrac{\pi}{2},\tfrac{\pi}{2},\pi)$.

4$f$ electron state of CeB6 is believed to be almost localized in Ce3+-ion from the several experiments of the magnetic and transport properties. More directly, the Fermi-surface (FS) has been observed in the de Haas-van Alphen (dHvA) experiments[26, 27], the angle resolved photoemission spectroscopy (ARPES)[28, 29, 30] and the high-resolution photoemission tomography[31]. They has indicated an ellipsoidal FS centered at X point in the Brillouin zone (BZ) which is almost the same as that of LaB6 with the 4$f^{0}$ state. Hence the 4$f$ state in CeB6 is localized and hardly participates in the formation of FS.

In such a localized $f$ electron picture, Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction[42, 43, 44] plays an important role for the multipole ordering where the intersite coupling between the multipole moments of $f$ electrons is mediated by the itinerant $c$ band electrons[45, 46]. The RKKY model of CeB6 was proposed by Ohkawa[32, 33] firstly, and later developed by Shiina et al.,[34, 35] where all 15 active multipole moments had been taken into account in correct symmetry, and reproduced the experimental $T$-$H$ phase diagram where only nearest neighbor couplings and the largest $\Gamma_{5g}$ quadrupoles couplings were assumed. This assumption was discussed from the symmetry of the RKKY couplings[36, 37], but there was no explicit calculation for the signs and values of the couplings, and also no discussion about the long-range property of the RKKY multipole coupling of CeB6. Later Sakurai et al.,[38, 39] studied the RKKY multipole couplings of CeB6 microscopically such as the effect of the $f^{0}$ and $f^{2}$ intermediate states, $c$ band number dependence and the ratio of the $f$-$c$ mixing elements described by the Slater-Koster (SK) parameters, but plausible ordering moment types and wavevectors could not be obtained.

As is often discussed in the RKKY mechanism, the $c$ band states and their couplings with the $f$ states in the realistic materials must be important for determining the ordering moment types and wavevectors. Therefore the microscopic description of the $c$ band states and $f$-$c$ mixing elements from the realistic bandstructure calculation is needed, though such studies are quite limited[47, 48]. In these studies[47, 48], the $c$ states is described by the Wannier orbitals obtained from the bandstructure calculation but the $f$-$c$ mixing elements using the calculation of the RKKY coupling are treated by the SK parameters only with the nearest neighbor sites, where several arbitrary parameters and assumptions are included. Hence more decisive and widely-applicable approach reflecting the individual material properties is highly desired.

In this paper, we study the electronic states of CeB6 and calculate the RKKY interaction based on the 74-orbital effective Wannier model derived from the bandstructure calculation directly. In Sec. 2, we calculate the bandstructures of CeB6 and LaB6 and construct the effective Wannier model of CeB6, and examine the quasi-particles states and their multipole fluctuations based on the renormalized Wannier model in Sec. 3. Next in Sec. 4, we formulate the present RKKY mechanism based on the multi-orbital Kondo lattice model with both of the $\Gamma_{8}$ quartet and 60 $c$ orbitals, and present the results of the RKKY multipole couplings for all moments as functions of the wavevector and intersite vectors. Finally we give the summary and discussion in Sec. 5.

2 Bandstructure calculation & Wannier model

2.1 DFT Bandstructure calculation of CeB6 & LaB6

First we calculate the electronic states of CeB6 and LaB6 by using the WIEN2k code[49, 50, 51], based on the framework of the density-functional theory (DFT) with the generalized gradient approximation (GGA)[52]. The SOC is fully included within the second variation approximation. The crystallographical parameters are the space group $Pm\bar{3}m$ (No. 221), the lattice constant $a=4.141{\rm\AA}$ and the internal coordinates $(x/a,y/b,z/c)=(0,0,0)$ for Ce and $(\frac{1}{2},\frac{1}{2},u)$ for B with $u=(\sqrt{2}-1)/2\sim 0.2071$ [53]. In self-consistent calculation, we use 156 $k$-points in the irreducible part of the simple cubic BZ, the muffin-tin radii $R_{\rm MT}=2.50~{}(1.62)$ a.u. for Ce (B) and the plane-wave cuttoff of $R_{\rm MT}K_{\rm max}=8$. For the calculation of LaB6, we use the same parameters of CeB6 but employ the GGA+$U$ method with $U=60$ eV for La-$f$ level so as to eliminate the $f$ weights in the $c$ bands, since we focus the pure $c$ band state of CeB6 not bulk property of LaB6.

The obtained bandstructures with the density-of-states (DOSs) and FSs are shown in Fig.1 for CeB6 [(a) & (b)] and LaB6 [(c) & (d)]. In CeB6, the large $f$ contribution due to the 14 $f$ spin-orbital states around Fermi energy ($E_{\rm F}$) is observed with the strong peak of DOS as shown in the right panel of Fig. 1 (a). On the other hands in LaB6 the $f$ states is absent in the bandstructure and DOS [Fig. 1 (b)] as expected due to the effect of the GGA+$U$. Except for the $f$ band states, the global bandstructures of CeB6 and LaB6 are closely resembled below and above $E_{\rm F}$. The calculated FSs of CeB6 and LaB6 are plotted in Figs. 1 (b) and (d), respectively. Three FSs are obtained from the 21st, 23rd and 25th bands for CeB6 while for LaB6 an ellipsoidal FS centered at X point slightly connected each other is obtained from the 21th-band. Here we note that all bands have two-folded degeneracy due to the time-reversal symmetry and two additional bands (1st and 2nd bands) are located in $E_{\rm F}-15$ eV (not shown) which are the lowest bands in the Wannier model in next subsection.

2.2 Construction of Wannier model for CeB6

Next we construct the 74-orbital effective Wannier model based on the maximally localized Wannier functions (MLWFs) method[54, 55, 56, 57, 58] from the DFT bandstructure of CeB6, where we prepare 14 $f$-states from Ce-$f$ (7 orbital $\times$ 2 spin) and 60 $c$-states from Ce-$d$ (5 orbital $\times$ 2 spin), Ce-$s$ (1 orbital $\times$ 2 spin), B-$p$ (6 site $\times$ 3 orbital $\times$ 2 spin) and B-$s$ (6 site $\times$ 1 orbital $\times$ 2 spin) as basis functions, and set considerably wide energy window in order to ensure the good localization of Wannier orbitals in the disentanglement procedure. The obtained bandstructure of the Wannier model is plotted in Figs. 1 (e) and (f) together with the DFT bandstructure of CeB6 (black), where the Wannier model is well reproduced the DFT bandstructure upto $E_{\rm F}+4$ eV and the shapes of the Wannier orbitals are similar to the atomic-orbitals significantly.

The obtained model can be written by the following tight-binding (TB) Hamiltonian as,

[TABLE]

where $f_{im}^{\dagger}~{}(c_{i\ell}^{\dagger})$ is a creation operator for a $f~{}(c)$ electron with unit-cell $i$ and 14 (60) spin-orbital states $m~{}(\ell)$. Here 14 $f$ states of $m$ are represented by the CEF eigenstates as $\Gamma_{8}$ quartet and $\Gamma_{7}$ doublet with the total angular momentum $J=5/2$, and $\Gamma_{6}$, $\Gamma_{7}$ doublets and $\Gamma_{8}$ quartet with $J=7/2$. The $f$-$f$ ($c$-$c$) matrix element of $h_{im,jm^{\prime}{}}^{ff}~{}(h_{i\ell,j\ell^{\prime}{}}^{cc})$ includes the $f~{}(c)$ energy levels, SOC couplings, CEF splittings and $f$-$f$ ($c$-$c$) hopping integrals, and $V_{im,j\ell}$ is the $f$-$c$ mixing element which is finite only for the intersite terms due to the inversion symmetry. The wavevector $\bm{k}$-representation of $H_{\rm TB}$ is given by,

[TABLE]

where $\varepsilon_{\bm{k}s}$ is the eigenenergy with $\bm{k}$ and band-index $s$ and $a_{\bm{k}s}^{\dagger}$ is a creation operator for a electron with $\bm{k},s$, which is transformed into $m$ and $\ell$ states as $\displaystyle a_{\bm{k}s}=\sum_{m}u_{\bm{k}sm}f_{\bm{k}m}+\sum_{\ell}u_{\bm{k}s\ell}c_{\bm{k}\ell}$ where $u_{\bm{k}sm}~{}(u_{\bm{k}s\ell})$ is the eigenvector component of $m~{}(\ell)$ state.

Several atomic parameters are obtained from the Wannier model, such as the SOC splitting for Ce-$4f$ between $J=5/2$ and $J=7/2$ states $\Delta_{\rm SOC}=0.33$ eV close to the experimental value of 3000 K, the atomic CEF splitting between $\Gamma_{8}$ and $\Gamma_{7}$, $\Delta_{\rm CEF}=8.2$ meV which is smaller than the experimental value of 540 K (=$46$ meV). The $f~{}(c)$ electron number per unit-cell is $n^{f}=1.24$ ($n^{c}=20.86$) and the total number is $n_{tot}=22$. All the $f$ electron number for each CEF state becomes finite where $n^{f}(\Gamma_{8})=0.634$ and $n^{f}(\Gamma_{7})=0.205$ for $J=5/2$ and $n^{f}(\Gamma_{6})=0.098$, $n^{f}(\Gamma_{7})=0.088$, and $n^{f}(\Gamma_{8})=0.216$ for $J=7/2$, due to the considerable $f$-$f$ hopping and $f$-$c$ mixing, which is indispensable within the DFT-based calculation.

3 Quasi-particle band states & Multipole fluctuations

3.1 Renormalized tight-binding model

As mentioned in Sec. 2, the $f$ electron state obtained here is fully itinerant and differs from the expected situation in the real material as $n^{f}(\Gamma_{8})\sim 1$. In this section, we examine the change of the electronic states and its multipole fluctuations from the itinerant $f$ band state to the localized $f$ state when $n^{f}(\Gamma_{8})=1$ in the realistic CeB6 bandstructures. For this purpose, we introduce a renormalization factor $Z_{m}^{f}$, which is explicitly derived from the Fermi-liquid (FL) theory [59], where the many-body correlation effect of the local $f$-$f$ Coulomb interaction is introduced through the self-energy $\Sigma_{m}^{f}(\bm{k},\varepsilon)$ which is almost local $\Sigma_{m}^{f}(\bm{k},\varepsilon)=\Sigma_{m}^{f}(\varepsilon)$ and can be expanded around $\varepsilon=0$ by the following form,

[TABLE]

where $Z_{m}^{f}$ corresponds to an inverse mass-enhancement $m/m^{*}$, and $\Delta\varepsilon_{m}^{f}$ and $\gamma_{m}$ are a shift of the $f$ energy-level and a dumping rate of the quasi-particles respectively. Hence in the itinerant quasi-particle picture, our original model of $H_{\rm TB}$ is renormalized by $Z_{m}^{f}$ and $\Delta\varepsilon_{m}$, yielding the renormalized tight-binding model $H_{\rm RTB}$ as explicitly given by,

[TABLE]

where the renormalized $f$-$f$ ($f$-$c$) matrix elements $\tilde{h}_{im,jm^{\prime}{}}^{ff}~{}(\tilde{V}_{im,j\ell})$ are written as,

[TABLE]

where $\varepsilon_{m}^{f}$ is a $f$ energy-level of the CEF state $m$, where $h^{ff}_{im,im^{\prime}{}}=0$ for $m\neq m^{\prime}{}$, and the $m$-dependence of $Z_{m}^{f}$ and $\Delta\varepsilon_{m}^{f}$ are dropped for simplicity as $Z_{m}^{f}=Z_{f}$ and $\Delta\varepsilon_{m}^{f}=\Delta\varepsilon_{f}$, where $\Delta\varepsilon_{f}$ is set to $\Delta\varepsilon^{f}=0.27$ eV so as to satisfy $n^{f}(\Gamma_{8})=1$ and $n^{c}=21$ at $Z_{f}=0$. Hence the $f$-$f$ ($f$-$c$) hopping elements are renormalized by $Z_{f}~{}(\sqrt{Z_{f}})$. Throughout the calculation, we determine a chemical potential $\mu$ so as to satisfy $n_{tot}=n^{f}+n^{c}=22$ with 643 $\bm{k}$-meshes in the entire BZ.

3.2 Renormalized electronic states

Figure 2 (a) shows the $Z_{f}$-dependence of $f$ electron number per CEF eigenstates $n^{f}(\Gamma)$ with $T=0.002$ eV, where $\Gamma=\Gamma_{8}^{J=5/2}$, $\Gamma_{7}^{J=5/2}$ $\Gamma_{6}^{J=7/2}$, $\Gamma_{7}^{J=7/2}$ and $\Gamma_{8}^{J=7/2}$, and $Z_{f}=1$ ($Z_{f}=0$) corresponds to the DFT-band (localized $f$) limit. With decreasing $Z_{f}$, $n^{f}(\Gamma_{8}^{J=5/2})$ increases and finally becomes $n^{f}(\Gamma_{8}^{J=5/2})=1$ when $Z_{f}=0$ while $n^{f}(\Gamma_{7}^{J=5/2})$ and all other $n^{f}(\Gamma^{J=7/2})$ decrease and reach zero at $Z_{f}=0$. The change of $n^{f}(\Gamma)$ is rapidly for $Z_{f}\stackrel{{\scriptstyle<}}{{{}_{\sim}}}0.1$ where the effective mass-enhancement reaches $m^{*}/m\stackrel{{\scriptstyle>}}{{{}_{\sim}}}10$. The $f$ electron magnetization $\displaystyle m_{\rm tot}=m_{\rm spin}+m_{\rm orb}$ as a function of $Z_{f}$ is also plotted in Fig. 2 (b) together with its spin, orbital and $J_{z}$-components $m_{\rm spin}$, $m_{\rm orb}$ and $m_{jz}$, respectively, where the magnetic field is applied along the $z$-direction with $h=\mu_{\rm B}H=0.004$ eV. The Zeemann Hamiltonian is given by $H_{Z}=\left(\sigma^{z}+\ell^{z}\right)h$, and $m_{\rm spin}$ and $m_{\rm orb}$ are explicitly written as,

[TABLE]

where $\sigma^{z}~{}(\ell^{z})$ is a $z$-component of spin Pauli (orbital angular momentum) matrix for $m$-basis and $f(x)$ is the Fermi distribution function $\displaystyle f(x)=\frac{1}{e^{\beta(x-\mu)}+1}$. With decreasing $Z_{f}$, $m_{\rm tot}$ increases and finally reaches the saturated value of the $\Gamma_{8}$ state as 1.5$\mu_{\rm B}$ together with an opposite sign between $m_{\rm orb}$ and $m_{\rm spin}$ due to the SOC effect.

The $T$-dependence of the magnetization $m_{\rm tot}$ and inverse magnetization $h/m_{\rm tot}$ for several values of $Z_{f}$ are plotted in Figs. 2 (c) and (d) respectively. For $Z_{f}=1\sim 0.3$ the weak $T$-dependence of $m_{\rm tot}$ is observed as a Pauli paramagnetic behavior of the itinerant $f$ electron, while for $Z_{f}\stackrel{{\scriptstyle<}}{{{}_{\sim}}}0.1$ $m_{\rm tot}$ increases with decreasing $T$, exhibiting the Curie paramagnetic behavior of the localized $f$ electron $m_{\rm tot}/h\sim 1/T$, which is more clearly observed in the inverse magnetization $h/m_{\rm tot}$ with a linear $T$-dependence. In such situations for $Z_{f}=0.1\sim 0.01$, the electronic state is similar to the purely localized $f$ electron state on a single Ce-ion usually analyzed in the experiments. However in this study the $f$-$c$ mixings are still finite and the quasi-particle hybridization bands are formed with the wide-bandwidth $c$ band dispersion having the ellipsoidal FS observed ARPES of CeB6.

Next we check such renormalized bandstructures for several values of $Z_{f}$ as shown in Figs. 3 (a)-(d) together with the DFT-bandstructure of CeB6 [Figs. 3 (a)-(c)] and the LaB6 GGA+$U$ band without $f$ weights [Fig. 3 (d)]. From $Z_{f}=0.5$ [Fig. 3 (a)] to $Z_{f}=0.1$ [Fig. 3 (b)], the whole bandstructures are still close to the DFT-band of CeB6 but their $f$ band-widths become narrow gradually, exhibiting a separation between the lower $J=5/2$ bands and higher $J=7/2$ bands. In Fig. 3 (d) with $Z_{f}=0.01$ corresponding to $m^{*}/m\sim 100$, the almost flatted $J=7/2$ bands, and $\Gamma_{7}$ and $\Gamma_{8}$ bands of $J=5/2$ are clearly observed around $E_{\rm F}$, and they slightly hybridize with the wide-bandwidth $c$ bands expanding from the X point in the BZ. Interestingly, the $c$ band dispersion with $Z_{f}=0.01$ (red) is almost overlapping the LaB6 band with the GGA+$U$ (black) as shown in Fig. 3 (d) except for the highly-flatted $f$ bands, resulting in the formation of almost the same FS of LaB6. Hence the $c$ bands of CeB6 with almost localized $f$ electron state coincides that of LaB6 without the La-$f$ contribution, and then their FS is also almost the same as that of LaB6 as shown in Fig. 1 (d). These results strongly support the localized $f$ electron picture for CeB6, and then the approach based on the periodic Anderson model and its perturbation w. r. t. the $f$-$c$ mixing is expected to giving a good starting point for treating this system.

3.3 Multipole fluctuations in the quasi-particle bands

Before going to the calculation of the RKKY interaction, we examine the multipole fluctuations under the renormalized $f$ bands on CeB6 by calculating the multipole susceptibility $\chi_{O_{\Gamma}}(\bm{q})$ with the multipole operator $O_{\Gamma}$ shown in Table 1 and the wavevector $\bm{q}$ which is given by,

[TABLE]

where $O_{\Gamma}=\sum_{mm^{\prime}{}}O_{mm^{\prime}{}}^{\Gamma}f_{im}^{\dagger}f_{im^{\prime}{}}$ and $O_{mm^{\prime}{}}^{\Gamma}$ is the normalized 4$\times$4 matrix element of $O_{\Gamma}$ in $\Gamma_{8}$ subspace[3], and $\chi_{m_{1}m_{2}m_{3}m_{4}}(\bm{q})$ is the irreducible $f$ electron susceptibility which depends on the distribution of $f$ states in the bandstructures through the renormalized $f$-$c$ mixing $\sqrt{Z_{f}}V_{im,j\ell}$.

Figure 4 (a) shows the $\bm{q}$-dependence of $\chi_{O_{{}_{\Gamma}}}(\bm{q})$ for each multipole moment with $Z_{f}=1$ and $\Delta\varepsilon_{f}=0$ eV corresponding to the DFT band limit as shown in Figs. 1 (e) and (f). The obtained $\bm{q}$-dependence is considerably weak and the explicit values of $\chi_{O_{{}_{\Gamma}}}(\bm{q})$ fall within the only small range $\chi_{O_{{}_{\Gamma}}}(\bm{q})=4\sim 5$ eV*-1* for all multipole moments and wavevectors $\bm{q}$. Among them the $\Gamma_{4u}^{(1)}$ magnetic multipole $\sigma^{z}$ susceptibility, where $\sigma^{x}$ and $\sigma^{y}$ are degenerate with $\sigma^{z}$, is barely large for the incommensurate wavevector around $\bm{q}=(\pi,\pi,\pi)$, while the $\Gamma_{5g}$ quadrupole $O_{xy}$ and the $\Gamma_{2u}$ ocutupole $T_{xyz}$ susceptibilities does not become large for the AFQ wavevector $\bm{q}=(\pi,\pi,\pi)$. The weak-$\bm{q}$ dependence of $\chi_{O_{{}_{\Gamma}}}(\bm{q})$ becomes more notable for the almost localized $f$ case with $Z_{f}=0.01$ and $\Delta\varepsilon_{f}=0.27$ eV as shown in Fig. 4 (b), where $\sigma^{z}$ becomes also maximum but its wavevector shifts to $\bm{q}=(0,0,0)$ as shown in the inset of Fig. 4 (b).

In such a situation, the actual value of $\chi_{O_{{}_{\Gamma}}}(\bm{q})$ becomes huge, where the extremely narrow $\Gamma_{8}$ bands are located in the very near and just above $E_{\rm F}$ with tiny $f$-$c$ mixing, and then the hybridized band $\varepsilon_{\bm{k}s}$ is highly degenerate for wide-range of the BZ, giving rise to the sizable enhancement of the Lindhard function of in Eq. (13). As far as such $\bm{q}$-independent $\chi_{O_{{}_{\Gamma}}}(\bm{q})$, it is difficult to describe the development of the $(\pi,\pi,\pi)$-AFQ mode with $(O_{yz},O_{zx},O_{xy})$ by the perturbation of the $f$-$f$ Coulomb interaction such as the random phase approximation (RPA) and its extensions.

4 RKKY Interaction of CeB6

4.1 Derivation of RKKY Hamiltonian

Here we consider the RKKY interaction between the multipole moments of $\Gamma_{8}$ quartet. For this purpose, we eliminate the $f$ energy-levels but use the $f$-$c$ mixing of the original Wannier model. The $c$ bandstructure for the calculation of RKKY couplings is shown in Figs. 5 (a) and (b), which is almost the same as that of LaB6 as mentioned in Sec. 3 and is compared to the strongly renormalized quasi-particles case with $Z_{f}=0.01$ as shown in Figs. 5 (c) and (d). During the calculation, $\mu$ is determined so as to keep $n_{tot}=n^{c}=21$ and $T$ is set to $T=0.005$ eV.

The multi-orbital Kondo lattice Hamiltonian for the present model is given by,

[TABLE]

where $m$ represents 4-states in $\Gamma_{8}$ quartet $\Ket{m}=\Ket{1}\sim\Ket{4}$ with an degenerate energy-level $\tilde{\varepsilon}_{m}^{f}$, which are given with the $J_{z}$-base of $J=5/2$ $\Ket{JM}$ explicitly as, $\Ket{1}=-\sqrt{\frac{1}{6}}\Ket{+\frac{3}{2}}-\sqrt{\frac{5}{6}}\Ket{-\frac{5}{2}}$, $\Ket{2}=\Ket{+\frac{1}{2}}$, $\Ket{3}=-\Ket{-\frac{1}{2}}$ and $\Ket{4}=\sqrt{\frac{1}{6}}\Ket{-\frac{3}{2}}+\sqrt{\frac{5}{6}}\Ket{+\frac{5}{2}}$. Here we note that $\tilde{\varepsilon}_{m}^{f}=\varepsilon_{m}^{f}+\Delta\varepsilon_{m}^{f}$ in Eq. (8) for $m=1\sim 4$ of $\Gamma_{8}$ and is pushed up from the bare $f$ energy-level $\varepsilon_{\Gamma_{8}}^{f}$ due to the DFT Hartree and GGA potentials which is of the order of a few eV. The $c$-$c$ matrix element $h_{\ell\ell^{\prime}{}}^{cc}(\bm{k})$ includes the $c$ orbital energy $\varepsilon_{\bm{k}\ell}^{c}$ for $\ell=\ell^{\prime}{}$ and the $c$-$c$ hopping $t_{\ell\ell^{\prime}{}}^{cc}(\bm{k})$ for $\ell\neq\ell^{\prime}{}$. The second term is rewritten by the $c$ band eigenstate $\displaystyle c_{\bm{k}s}=\sum_{\ell}u_{\bm{k}s\ell}^{c}c_{\bm{k}\ell}$ with the eigenenergy $\varepsilon_{\bm{k}s}$ and eigenvector $u_{\bm{k}s\ell}^{c}$.

The Kondo coupling $J_{imm^{\prime}{}}^{\bm{k}\ell,\bm{k}^{\prime}{}\ell^{\prime}{}}$ in the third term consists of the $f^{0}$- and $f^{2}$-intermediate process. In this paper, we take simple two assumptions for $J_{imm^{\prime}{}}^{\bm{k}\ell,\bm{k}^{\prime}{}\ell^{\prime}{}}$; (1) only $f^{0}$-process is considered and the contribution of $f^{2}$-process is same as that of $f^{0}$-process and, (2) the scattered $c$ orbital energies are fixed to $\mu$, namely $\varepsilon_{\bm{k}\ell}^{c},~{}\varepsilon_{\bm{k}^{\prime}{}\ell^{\prime}{}}^{c}\longrightarrow\mu$. Then the Kondo coupling $J_{imm^{\prime}{}}^{\bm{k}\ell,\bm{k}^{\prime}{}\ell^{\prime}{}}$ can be written by the following simple form,

[TABLE]

where the prefactor 2 comes from the assumption (1) and $V_{\bm{k}m\ell}$ is the $\bm{k}$-represented $f$-$c$ mixing element in Eq. (2).

The RKKY Hamiltonian can be obtained from the second-order perturbation w. r. t. the third term of $H_{\rm MKL}$ together with the thermal average for the $c$ states. The final form is given by,

[TABLE]

where $K_{m_{1}m_{2}m_{3}m_{4}}(\bm{R}_{ij})$ is the RKKY coupling between the states $\{m_{1},m_{2}\}$ at the unit-cell $\bm{R}_{i}$ and the states $\{m_{3},m_{4}\}$ at $\bm{R}_{j}$ and $\Braket{ij}$ represents a summation for the intercell vectors $\bm{R}_{ij}=\bm{R}_{i}-\bm{R}_{j}$. The key quantity $K_{m_{1}m_{2}m_{3}m_{4}}(\bm{q})$ is given by,

[TABLE]

which consists of a square of the energy denominator $(\mu-\varepsilon_{\Gamma_{8}}^{f})^{-2}$, 4-producted $f$-$c$ mixings and the $c$ band eigenvectors, and the Lindhard function with $\varepsilon_{\bm{k}s}$. Thus it has $4^{4}=256$ components of $f$-basis $\{m_{1},m_{2},m_{3},m_{4}\}$ for each $\bm{q}$, and has to be summed for the $c$ orbitals $\{\ell_{1},\ell_{2},\ell_{3},\ell_{4}\}$ (604) and the band-indexes $\{s,s^{\prime}{}\}$ ($60^{2}$). Then we introduce a $f$-$c$ mixing matrix $v_{\bm{k}s}^{mm^{\prime}{}}$ between $\{m,m^{\prime}{}\}$ via the $c$ band state with $\bm{k},s$ as follows,

[TABLE]

which includes whole information about the $f$ state scattering between $\{m,m^{\prime}{}\}$ through the $c$ state with $\bm{k},s$, and has only $4^{2}=16$ components of $\{m,m^{\prime}{}\}$ for each $\bm{k},s$ with a summation for $\{\ell,\ell^{\prime}{}\}$ (602). Hence once we calculate $v_{\bm{k}s}^{mm^{\prime}{}}$, $K_{m_{1}m_{2}m_{3}m_{4}}(\bm{q})$ can be easily obtained by the following compact form,

[TABLE]

This expression helps us calculate all the contributions of the 60 $c$ electron charge and/or orbital fluctuations to the RKKY multipole couplings.

In order to search the actual multipole ordering, we employ the mean-field (MF) approximation w. r. t. the multipole operator $O_{\Gamma}$, resulting in the MF Hamiltonian as follows,

[TABLE]

where the multipole coupling $\overline{K}_{O_{\Gamma}}(\bm{q})=K_{O_{\Gamma}}(\bm{q})-K_{O_{\Gamma}}^{\rm loc}$ and $K_{O_{\Gamma}}(\bm{q})$ and the MF order parameter $\overline{O}_{\Gamma}(\bm{q})$ are given by,

[TABLE]

where $K_{O_{\Gamma}}^{\rm loc}=(1/N)\sum_{\bm{q}}K_{O_{\Gamma}}(\bm{q})$ and $\overline{O}_{\Gamma}(\bm{R}_{i})$ is MF multipole order parameter defined at the unit-cell vector $\bm{R}_{i}$. Then the MF multipole susceptibility $\chi_{O_{\Gamma}}^{\rm MF}(\bm{q})$ is written by,

[TABLE]

which is enhanced towards the multipole ordering instability for the ordering moment $O_{\Gamma}$ and wavevector $\bm{q}$, and finally diverges at a critical point of the multipole ordering transition temperature $T=T_{O_{\Gamma}}$ where $\chi_{O_{\Gamma}}(\bm{q})\overline{K}_{O_{\Gamma}}(\bm{q})$ reaches unity. The $\bm{q}$-dependence of $\chi_{O_{\Gamma}}(\bm{q})$ is weak as shown in Sec. 3, and then the sign and maximum value of $\overline{K}_{O_{\Gamma}}(\bm{q})$ determines the multipole ordering moment and wavevector for any given $T$. Hereafter we set $\mu-\varepsilon_{\Gamma_{8}}^{f}=1$ eV for simplicity, since this factor is independent of $\bm{q}$ and $m$, and hence does no affect the ordering type and wavevector, whose effect is discussed in Sec. 5.

4.2 $\bm{q}$-dependence of RKKY coupling $\overline{K}_{O_{\Gamma}}(\bm{q})$

The obtained RKKY multipole couplings $\overline{K}_{O_{\Gamma}}(\bm{q})$ for several multipole moments along the high symmetry line in the BZ are plotted as shown in Figs. 6 (a)-(d), where the positive (negative) coupling for a certain multipole $O_{\Gamma}$ and wavevector $\bm{q}$ enhances (suppresses) the corresponding multipole fluctuation as explained in Eq. (24), and its positive maximum value gives a leading multipole ordering mode. The obtained results for the leading multipole ordering modes upto the 10th largest coupling are summarized in Table 2.

The couplings of the $\Gamma_{5g}$ quadrupoles $(O_{yz},O_{zx},O_{xy})$ for $\bm{q}=(\pi,\pi,\pi)$ become largest among all moments and $\bm{q}$, which perfectly corresponds to the AFQ ordering of CeB6. In addition, $\Gamma_{2u}$ ocutupole $T_{xyz}$ coupling is quite large and comparable to the $\Gamma_{5g}$ quadrupoles with the same wavevector as shown in Fig. 6 (a) but slightly small within the present calculation accuracy as shown in Table 2, which seems to be the same value from the previous discussions[36, 37] where $O_{xy}$ and $T_{xyz}$ have almost same matrix elements and yield the similar fluctuations in phase I. Furthermore the quadupoles $(O_{yz},O_{zx},O_{xy})$ and octupole $T_{xyz}$ couplings also take a substantial peak for $\bm{q}=(0,0,0)$ as shown in Figs. 6 (a) and (d) and correspond to the elastic softening of $C_{44}$[15, 16, 17, 18].

The next largest coupling is the $\Gamma_{5u}$ octupole $\zeta^{x}$ at $\bm{q}=(\pi,0,0)$ [Figs. 6 (a) & (d)] which is degenerate for $\zeta^{y}~{}[\zeta^{z}]$ octupole at $\bm{q}=(0,\pi,0)~{}[(0,0,\pi)]$ due to the cubic symmetry. The role of the $\Gamma_{5u}$ octupoles $(\zeta^{x},\zeta^{y},\zeta^{z})$ is also discussed for the phase IV observed in the La-doping system CexLa1-xB6 with $x<0.8$ [60, 61, 62, 63, 64, 65, 66], where the $\bm{q}=(\pi,\pi,\pi)$ antiferro-octupolar (AFO) ordering of $(\zeta^{x}+\zeta^{y}+\zeta^{z})/\sqrt{3}$ is considered to be a possible mode. In contrast, the present theory suggests the $\Gamma_{5u}$ AFO with the domained structure of $\zeta^{x}$, $\zeta^{y}$ and $\zeta^{z}$ for $\bm{q}=(\pi,0,0)$, $(0,\pi,0)$ and $(0,0,\pi)$, respectively, and this point will be discussed in the next subsection.

In addition to this, the $\Gamma_{3g}$ quadrupole $O_{v}=O_{x^{2}-y^{2}}$ coupling is quite large for $\bm{q}=(0,0,\pi)$ (not shown) and becomes similar value of the octupole coupling $\zeta^{z}$ as shown in Table 2, which is also degenerate for the rotated moments to the each principle-axis $O_{y^{2}-z^{2}}$ and $O_{z^{2}-x^{2}}$. This is namely the $\Gamma_{3g}$-AFQ mode where the moment directions and wavevectors are perpendicular such as the multipole moments of $O_{y^{2}-z^{2}}$, $O_{z^{2}-x^{2}}$ and $O_{x^{2}-y^{2}}$ with the corresponding wavevectors for $\bm{q}=(\pi,0,0)$, $(0,\pi,0)$ and $(0,0,\pi)$ respectively.

The $\Gamma_{4u}$ magnetic multipole couplings of $(\sigma^{x},\sigma^{y},\sigma^{z})$ and $(\eta^{x},\eta^{y},\eta^{z})$ are plotted in Fig. 4 (d) and their maximum values in $\bm{q}$-space are smaller than that of the quadrupole and octupole couplings as shown in Table 2, where the $\Gamma_{4u}$ maximum peak values are less than half of the first leading peak value of the $\Gamma_{5g}$-$(\pi,\pi,\pi)$. At the AFM ordering vectors for phase III, $\bm{Q}_{1}$ and $\bm{Q}_{2}$, the couplings of the magnetic multipoles $(\sigma^{x},\sigma^{y},\sigma^{z})$ have small peaks as shown in Fig. 4 (e) and they shall be enhanced and dominant only when the system enters into phase II, which is not discussed in the present paper.

As usually discussed in the itinerant $f$ electron picture with the multi-orbital Hubbard model[67, 68, 69], the weak coupling theory like the RPA and its extensions yields largely enhanced magnetic multipole (spin) fluctuations which become always larger than the nonmagnetic multipole (orbital) fluctuations like $\Gamma_{5g}$ multipoles here. As for CeB6, the $f$ electron itself is already localized at each Ce site and the remained magnetic and nonmagnetic multipole moments interact with the RKKY intersite couplings, where the magnetic multipole coupling does not necessarily dominate over the nonmagnetic one, since the dominant RKKY coupling is determined by the detail of the $f$-$c$ mixing and the meditating $c$ electron charge and/or orbital fluctuations.

Here we note the $c$ electron charge and orbital fluctuations and their contribution to the coupling $\overline{K}_{O_{\Gamma}}(\bm{q})$. By changing summation for the $c$ orbital-set in Eq. (18) and the band-index in Eq. (20), we have obtained that both effects of the Ce-$d_{eg}$ orbitals of $d_{x^{2}-y^{2}}$ and $d_{3z^{2}-r^{2}}$ distributed in the 21st and 23rd bands and the charge fluctuation of B6-molecule having a maximum at $\bm{q}=(\pi,\pi,\pi)$ and large values along R-M line play significant roles for the $\Gamma_{5g}$ AFQ mode. In particular, we observe a non-negligible contribution from the 23rd band which does not have FS but is very close to $E_{\rm F}$ along the $\Gamma$-M direction as shown in Fig. 5 (b). The explicit results and further analysis of such $c$ electron contributions to the multipole couplings will be presented in elsewhere.

4.3 $\bm{R}_{ij}$-dependence of RKKY coupling $K_{O_{\Gamma}}(\bm{R}_{ij})$

In general, the RKKY interaction is known to have long-range and oscillating features discussed in the early studies[45, 46]. For the multipole ordering of CeB6, however, the coupling is limited only in the nearest neighbor terms in the previous studies[32, 33, 34, 35, 36, 37, 38, 39]. In contrast, the present formalism provides the real space dependent couplings $K_{O_{\Gamma}}(\bm{R}_{ij})$ which is explicitly written as

[TABLE]

where $K_{m_{1}m_{2}m_{3}m_{4}}(\bm{R}_{ij})$ is given in Eq. (17).

Figure 7 shows the site-dependence of the RKKY multipole couplings $K_{O_{\Gamma}}(\bm{R}_{ij})$ with the intersite vector $\bm{R}_{ij}$ upto 20-th neighbor sites as shown in Table 7, where the positive (negative) sign corresponds to the ferro (antiffero) coupling for each neighboring site.

As shown in Fig. 7 (a), the $\Gamma_{4u}^{(1)}$ magnetic multipole $(\sigma^{x},\sigma^{y},\sigma^{z})$ couplings exhibit several sign changes with a few site-intervals and degeneracy due to the symmetry of paramagnetic phase for all $n$-th neighbors, where $\sigma^{x}$, $\sigma^{y}$ and $\sigma^{z}$ corresponds to the $x,~{}y$- and $z$-moment direction respectively. We also confirm that a monopole operator $I$ defined as a unit matrix for $\Gamma_{8}$-basis is also degenerate with $(\sigma^{x},\sigma^{y},\sigma^{z})$ possessing the same oscillating feature. The couplings of $I$ and $(\sigma^{x},\sigma^{y},\sigma^{z})$ are isotropic, where for example they have the same value for 6 first neighbor sites $\bm{R}_{ij}=(\pm a,0,0)$, $(0,\pm a,0)$ and $(0,\pm a,0)$.

The $\Gamma_{5g}$ and $\Gamma_{2u}$ multipoles couplings which gives the leading AFQ mode show staggered and isotropic behaviors, where the first, second and third neighbor couplings show positive, negative and positive signs respectively as shown in Figs. 7 (b) [$T_{xyz}$] and (c) [$(O_{yz},O_{zx},O_{xy})$], which clearly enhance cooperatively the antiferro ordering. In particular, the second neighbor coupling becomes largest with positive sign, which also enhances the AFQ mode as a main driving force.

On the other hands, the main origin of the $(\pi,0,0)$-AFO with $\Gamma_{5u}$ octupole $\zeta^{x}$ is the anisotropic first neighbor couplings as shown in Fig. 5 (e), where the coupling of which the intersite vector and moment direction are parallel (perpendicular) each other has negative (positive) sign such as $K^{x}_{{\rm 5}u}(\bm{R}_{x})<0$ and $K^{y}_{{\rm 5}u}(\bm{R}_{x}),~{}K^{z}_{{\rm 5}u}(\bm{R}_{x})>0$ for $\bm{R}_{x}=(a,0,0)$, which yields $K^{x}_{{\rm 5}u}(\bm{R}_{y}),~{}K^{x}_{{\rm 5}u}(\bm{R}_{z})>0$ for the perpendicular first neighbors $\bm{R}_{y}=(0,a,0)$ and $\bm{R}_{z}=(0,0,a)$, resulting in the enhancement of the $\bm{q}=(\pi,0,0)$ mode by $K^{x}_{{\rm 5}u}(\bm{q})=2K^{x}_{{\rm 5}u}(\bm{R}_{x})\cos q_{x}+2K^{x}_{{\rm 5}u}(\bm{R}_{y})\cos q_{y}+2K^{x}_{{\rm 5}u}(\bm{R}_{z})\cos q_{z}$.

With this matter in mind, we might be able to explain the doping phase diagram of CexLa1-xB6, where the doping effect is simply treated by the reduction of the multipole coupling depending the coordination number for each site, since a La-substituted site has no multipole moment. Consequently, the first leading $(\pi,\pi,\pi)$-AFQ mode in $x=1$ decreases with decreasing $x$ more rapidly than the second $(\pi,0,0)$-AFO mode due to the difference of the coordination numbers : 6 first neighbors and 12 second neighbors for $(\pi,\pi,\pi)$, while 2 parallel first neighbors and 4 perpendicular first neighbors and so on for $(\pi,0,0)$. This turnover of the dominant mode may be consistent with the recent inelastic neutron scattering (INS) experiments in the La-doping system[25] where the $\bm{q}=(\pi,0,0)$ intensity is developed and becomes dominant mode for $x<0.8$.

As well as the $\Gamma_{5u}$ couplings, the $\Gamma_{3g}$ quadrupoles $(O_{u},O_{v})$ and $\Gamma_{4u}^{(2)}$ multipole $(\eta^{x},\eta^{y},\eta^{z})$ couplings exhibit the anisotropic behavior as shown in Figs. 7 (b) and (d). The obtained results of such long-range, oscillating, isotropic or anisotropic behaviors depending on the multipole moments seem to be worthwhile to study from now on.

5 Summary and Discussion

In summary, we study the electronic states of CeB6 and perform a direct calculation of the RKKY interaction based on the 74-orbital effective Wannier model derived from the bandstructure calculation and obtain the following results.

(1) When the $f$-$f$ hopping and $f$-$c$ mixing of the Wannier model are suppressed by the renormalization factor $Z_{f}$ based on the FL theory, the quasi-particle band states are observed where the fully dispersionless $f$ bands slightly hybridize with the wide-bandwidth $c$ bands, which is almost the same as the GGA+$U$-band of LaB6 having a single ellipsoidal FS centered at X point. This is in good agreement with the recent ARPES[28, 29, 30, 31] and early dHvA experiments[26, 27] and strongly supports the localized $f$ electron picture for CeB6.

(2) By using the LaB6-like $c$ band states together with the $f$-$c$ mixing elements of the CeB6 Wannier model, we calculate the RKKY couplings for all active multipole moments in $\Gamma_{8}$ subspace explicitly as functions of wavevector $\bm{q}$ and inter-cell vector $\bm{R}_{ij}$, where we derive a useful expression in order to treat all 60 $c$-orbital contributions.

(3) The couplings of the $\Gamma_{5g}$ quadrupole $(O_{yz},O_{zx},O_{xy})$ together with the $\Gamma_{2u}$ octupole $T_{xyz}$ are highly enhanced for $\bm{q}=(\pi,\pi,\pi)$ as the 1st and 2nd leading modes, and $\bm{q}=(0,0,0)$ as the 5th and 6th leading modes, where the former explains the AFQ ordering of the phase II and the latter corresponds to the elastic softening of $C_{44}$. The 3rd (4th) leading mode is the $\Gamma_{5u}$-AFO ($\Gamma_{3g}$-AFQ) of $\zeta^{x}$, $\zeta^{y}$ and $\zeta^{z}$ ($O_{y^{2}-z^{2}}$, $O_{z^{2}-x^{2}}$ and $O_{x^{2}-y^{2}}$) with the corresponding wavevectors for $\bm{q}=(\pi,0,0)$, $(0,\pi,0)$ and $(0,0,\pi)$ respectively, which are almost degenerate each other and differ from the discussed AFO-mode with $(\zeta^{x}+\zeta^{y}+\zeta^{z})/\sqrt{3}$ at $\bm{q}=(\pi,\pi,\pi)$[60, 61, 62, 63, 64, 65, 66].

(4) All the obtained RKKY couplings have long-range and oscillating behavior as a function of $\bm{R}_{ij}$, where the $\Gamma_{5g}$ quadrupole $(O_{yz},O_{zx},O_{xy})$ and $\Gamma_{2u}$ octupole $T_{xyz}$ couplings indicate the sign-reversing for each neighboring site and have a positive largest value at the second neighbor which cooperatively enhances the AFQ with $\bm{q}=(\pi,\pi,\pi)$, while for the second leading $\Gamma_{5u}$ AFO mode, the anisotropic first neighbor couplings are significant. This induces the leading mode shift with increasing the La-substitution rate $x$ in CexLa1-xB6 from the $(\pi,\pi,\pi)$-AFQ with $O_{xy}$ (phase II) to the $(\pi,0,0)$-AFO with $\zeta^{x}$ (phase IV) which may be also consistent with the $(\pi,0,0)$ peak in the INS data[25].

(5) The present approach can determine the possible type of the multipole moment and the ordering vector $\bm{q}$ definitely once the $c$ band states and $f$-$c$ mixings are given by the bandstructure calculation, which enables us to discuss the inherent feature and the concrete situation of actual compounds such as the changes of FSs, carrier densities, lattice constants and internal coordinates of atoms.

In this study, we take only the $f^{0}$-process and assume that the contribution from $f^{2}$-process is the same as that of $f^{0}$, since the $f^{0}$-process is fully one-body effect and directly obtained from the DFT-bandstructure calculation. As mentioned in Sec. 4, we take $\mu-\varepsilon_{\Gamma_{8}}^{f}=1$ eV, but the excitation energy from the $f^{1}$-stable to $f^{0}$-intermediate states in CeB6 is roughly estimated by $\mu-\varepsilon_{\Gamma_{8}}^{f}=2\sim 4$ eV[29, 31], so that our results of $\overline{K}_{O_{\Gamma}}$ obtained in Sec. 4 should be multiplied by a single reduction factor $\frac{1}{4}\sim\frac{1}{16}$, which yields the same order of the actual transition temperature of CeB6 as a few K, for example, the inter-quadrupole coupling value $K_{\Gamma_{5g}}^{\prime}{}=2.1$ K[17].

The explicit determination of the couplings and the transition temperatures needs the many-body energy difference between the ground and intermediate $f^{0}$ and $f^{2}$ states, where to what extent the many-body effect from the $f^{2}$ and more multiple $f$ processes changes the present result is an important question elucidating the multipole ordering system with different valence materials such as PrB6 and NdB6. The explicit calculation of the coupling including the $f^{2}$-process and/or more many-body contribution, and the whole phase diagram in $T$-$H$ plane will be presented in the subsequent paper.

As a complementary approach to the present localized $f$ electron treatment, the dynamical mean field theory[70] enabling to take account of the full local correlation effect and its extensions[71] including the intersite correlation could be valid for directly describing the fully localized $f$ states starting from the itinerant $f$ states and their multiple ordering phenomena including superconductivity[72]. The application of such many-body theory to the realistic materials and their comparison with present theory are also the essential future problems.

{acknowledgment}

We would like to thank Y. Ōno and Y. Iizuka for valuable comments and discussions.