#  RKKY Interactions and Multipole Order in Ab initio Wannier Model of CeCoSi

Takemi Yamada E-mail address<EMAIL_ADDRESS>Yuki Yanagi and Keisuke
Mitsumoto  Liberal Arts and Sciences Liberal Arts and Sciences Toyama
Prefectural University Toyama Prefectural University Imizu Imizu Toyama
939-0398 Toyama 939-0398 Japan Japan

###### Abstract

We calculate the RKKY interactions derived from ab initio calculations for the
intermetallic compound CeCoSi exhibiting the hidden nonmagnetic order at
$T_{0}$ and examine the instability towards possible multipole orders within
the random phase approximation. All 36 multipole interactions up to rank 5 are
investigated, and the maximum susceptibility exhibits an antiferro order with
$\bm{q}=\bm{0}$ for nonmagnetic multipoles of monopole $I$ and hexadecapole
$H_{0}$, yielding a charge imbalance of $f$ electrons at two Ce atoms in the
unit cell. The obtained order can explain some experiments.

The tetragonal intermetallic compound CeCoSi exhibits the so-called hidden
order (HO) at slightly higher temperature $T_{0}\simeq 12$ K[1, 2] than the
antiferromagnetic (AFM) order $T_{\rm N}=9.4$ K. The absence of magnetic order
with temperatures $T_{\rm N}<T<T_{0}$[3] has attracted much attention in terms
of non-magnetic higher-order multipole orderings, including cluster multipole
orderings[4, 5]. In the normal phase $T>T_{0}$, 4$f$ electrons show the
localized behavior in bulk experiments[2, 6], and the crystalline electric
field (CEF) level splits into three Kramers doublets (KDs) with two CEF
excitation energies around 100-150 K determined by the inelastic neutron
scattering[7]. The fact that the nonmagnetic HO occurs first is interesting,
since at very low temperatures only the lowest KD is active, and the onset of
magnetic order is generally expected. It is also noteworthy that the pressure
$P$-dependence of $T_{0}$ shows a dome-shaped pressure phase diagram similar
to the so-called Doniach phase diagram, increasing to about 40 K at $P\simeq
1.5$ GPa and disappearing at $P\simeq 2.2$ GPa[2], reminding us of the quantum
critical point of the nonmagnetic multipole.

Recent X-ray diffraction experiment have revealed that the electronic state at
$T<T_{0}$ is accompanied by a structural transition to a triclinic
structure[8], and the interpretation of $T_{0}$ as a ferro-quadrupole order of
$O_{zx}+O_{yz}$ is also supported by NMR measurement[9]. The magnetic phase
diagram[6, 10] at low temperatures including structural transitions is very
complicated depending on the direction of the magnetic field, and is still
being studied intensively toward a unified understanding. However, the
proposed quadrupole mechanisms[4, 5, 3, 10] require the quadrupole
interactions to be large enough to overcome the lowest KD fluctuations that
increases at low temperatures, but the magnitude and sign of the interactions
and their type dependence are only assumed phenomenologically. Therefore, the
microscopic analysis of multipole interactions reflecting the density
functional theory (DFT)-based bandstructure, fully including the conduction
($c$) bands and the mixing between localized $f$ and itinerant $c$ electrons
($c$-$f$ mixing), is crucial for the theoretical clarification of HO.

In this letter, we report a study of the Ruderman-Kittel-Kasuya-Yosida (RKKY)
interactions based on the multi-orbital Wannier model derived from the ab
initio calculation on CeCoSi, and the possible multipole order within the
random phase approximation (RPA). This approach is an extension of our work on
the typical localized multipole system CeB6[11, 12] and should be useful for
understanding CeCoSi.

Figure 1: (Color online) (a),(b) The $c$ bands of TB model (red) and the DFT
bands of LaCoSi (black) in (a) wide and (b) narrow energy ranges, where
horizontal green line is Fermi-level and the high-symmetry points are
$\Gamma~{}(0,0,0)$, X $(\pi,0,0)$, M $(\pi,\pi,0)$, Z $(0,0,\pi)$, R
$(\pi,0,\pi)$ and A $(\pi,\pi,\pi)$. (c),(d) The corresponding FSs for (c)
LaCoSi and (d) CeCoSi depicted by Fermisurfer[13].

First we perform band calculation for CeCoSi with the experimental lattice
parameters by using the DFT-based first-principles package WIEN2k[14], and
then construct an 84-orbital tight-binding (TB) model from the Wannier 90[15]
consisting of Ce-$4f,5d$, Co-$3d$ and Si-$3p,3s$ orbitals together with the
sub-lattice due to 2CeCoSi/unit-cell (uc) and spin degrees of freedom. The
computational details are presented in §1 and §2 of Supplemental Material
(SM)[16].

The obtained TB Hamiltonian is given by the following,

$\displaystyle\mathscr{H}_{\rm
TB}=\mathscr{H}_{c\textrm{-}c}+\mathscr{H}_{f\textrm{-}f}+\mathscr{H}_{c\textrm{-}f}$
(1)

where $\mathscr{H}_{c\textrm{-}c}~{}(\mathscr{H}_{f\textrm{-}f})$ is the
$c$-$c$ ($f$-$f$) term for $c~{}(f)$ electrons with 56 (28) orbitals and
$\mathscr{H}_{c\textrm{-}f}$ is the $c\textrm{-}f$ mixing term. By
diagonalizing $\mathscr{H}_{c\textrm{-}c}$ in the wavevector $\bm{k}$ space,
we can obtain the $c$ band dispersion $\varepsilon_{n\bm{k}}$ with band-index
$n$.

Figures 1(a) and (b) show $\varepsilon_{n\bm{k}}-\mu$ of the TB model for
CeCoSi together with the DFT bandstructure for LaCoSi in wide [narrow] energy
range as shown in Fig. 1(a) [Fig. 1(b)], where the chemical potential $\mu$ is
determined to satisfy the $c$ electron number $n_{c}=32$. The two bands are
similar especially for the occupied states below Fermi energy
($\varepsilon_{n\bm{k}}=\mu$), which indicates that
$\mathscr{H}_{c\textrm{-}c}$ can well describe the $c$ electron state of
LaCoSi almost without $f$ electrons and also is consistent with the angle-
resolved photoemission spectroscopy (ARPES) experiment on CeCoSi, where the
DFT bands of LaCoSi suggests relatively good correspondence with the
ARPES[17].

Table 1:  Classification of 36 multipole operators by IRR in the point group
$C_{4v}$ and TR symmetry [18, 19, 20], where even (odd) rank multipoles of
monopole (dipole), quadrupole (octupole) and hexadecapole (triakontadipole)
are denoted by $I~{}(J)$, $O~{}(T)$, and $H~{}(D)$, respectively[20] (see §3
of SM[16]).

IRR | TR-even | TR-odd
---|---|---
$\Gamma_{1}~{}(A_{1})$ | $I,~{}O_{u},~{}H_{0},~{}H_{4}$ |  $D_{4}$
$\Gamma_{2}~{}(A_{2})$ | $H_{z}^{\alpha}$ |  $J_{z},~{}T_{z}^{\alpha},~{}D_{z}^{1\alpha},D_{z}^{2\alpha}$
$\Gamma_{3}~{}(B_{1})$ | $O_{v},~{}H_{2}$ |  $T_{xyz},~{}D_{2}$
$\Gamma_{4}~{}(B_{2})$ | $O_{xy},~{}H_{z}^{\beta}$ |  $T_{z}^{\beta},~{}D_{z}^{\beta}$
$\Gamma_{5}~{}(E)$ |  |  $O_{zx},O_{yz}$
---
$H_{x}^{\alpha},~{}H_{y}^{\alpha},~{}H_{x}^{\beta},~{}H_{y}^{\beta}$
|
$J_{x},J_{y},~{}T_{x}^{\alpha},T_{y}^{\alpha},~{}T_{x}^{\beta},T_{y}^{\beta}$
---
$D_{x}^{1\alpha},~{}D_{y}^{1\alpha},~{}D_{x}^{2\alpha},~{}D_{y}^{2\alpha},~{}D_{x}^{\beta},~{}D_{y}^{\beta}$

On the other hand, there also are differences between the two bands in the
vicinity of $\varepsilon_{n\bm{k}}=\mu$ as seen in Fig. 1(b) and the
corresponding Fermi surfaces (FSs) in the tetragonal Brillouin zone (BZ) for
LaCoSi and CeCoSi are shown in Figs. 1(c) and (d), respectively. Three
independent FSs are obtained for the band-index 31st, 33rd and 35th, where all
bands have two-folded degeneracy due to the time-reversal (TR) symmetry. The
remarkable differences between the two bands can be seen especially in the
31st and 33rd FSs, where the 33rd FS forming elongated arms in the $\Gamma$-M
direction in LaCoSi [Fig. 1(c)] yields a square hole pocket centered at
$\Gamma$ as the 31st FS in CeCoSi with the separation of the $\Gamma$-M arms
[Fig. 1(d)]. The reason for these differences is that the $c$ bands of the TB
model has no complete $f$-electron component, whereas the DFT bands of LaCoSi
inevitably contains the $f$-component due to the finite $c$-$f$ mixing. For
the calculation of RKKY interactions as seen later, it is reasonable to
utilize the present $c$ states corresponding to the localized $f$-electron
limit without the $c$-$f$ mixing.

Next, we have extracted the one-body parameters for $f$ electrons from
$\mathscr{H}_{f\textrm{-}f}$, such as the spin-orbit coupling (SOC) constant
for $\lambda_{\rm SOC}=94.2$ meV close to the typical value of Ce atom. The
three KDs are obtained as $\Gamma_{6\pm}$ (0 meV), $\Gamma_{7a\pm}$ (9 meV)
and $\Gamma_{7b\pm}$ (23 meV), respectively, which disagree the experimental
CEF levels of $\Gamma_{7a\pm}$ (0 meV), $\Gamma_{7b\pm}$ (10 meV) and
$\Gamma_{6\pm}$ (14 meV)[7]. This contrasts with the results for CeB6, where
the ordering of the CEF levels in the experiment has been also obtained in the
DFT model[11, 12], possibly because the CEF splitting in the present system is
smaller than in CeB6. Here and hereafter, we adopt the experimental CEF levels
in the following, since the $c$-$f$ mixing being crucial for the RKKY
interaction is a one-body quantity between $f$ and $c$ electrons and should be
almost unaffected by the CEF levels.

Here we investigate the single-site multipole susceptibility for the multipole
$O_{\Gamma}$,

$\displaystyle\chi_{O_{{}_{\Gamma}}}^{(0)}=\sum_{m_{1}m_{2}m_{3}m_{4}}(O_{\Gamma})_{m_{1}m_{2}}\chi_{m_{1}m_{2}m_{3}m_{4}}^{(0)}(O_{\Gamma}^{\dagger})_{m_{3}m_{4}},$
(2)

where $\chi_{m_{1}m_{2}m_{3}m_{4}}^{(0)}$ is the single-site susceptibility
for 6 CEF states $\\{m_{i}\\}$,
$|1\pm\rangle=\mp\sqrt{1-w^{2}}|\pm\frac{5}{2}\rangle\pm
w|\mp\frac{3}{2}\rangle~{}(\Gamma_{7a\pm})$,
$|2\pm\rangle=w|\pm\frac{5}{2}\rangle+\sqrt{1-w^{2}}|\mp\frac{3}{2}\rangle~{}(\Gamma_{7b\pm})$,
$|3\pm\rangle=|\pm\frac{1}{2}\rangle~{}(\Gamma_{6\pm})$ with $w=0.95$[7]. 6
CEF states have $6\times 6=36$ multipoles up to rank 5 classified into the
irreducible representation (IRR) of the Ce site under the point group $C_{4v}$
as shown in Table 1.

Figure 2 shows the susceptibilities $\chi_{O_{{}_{\Gamma}}}^{(0)}$ for (a)
monopole and quadrupoles, (b) hexadecapoles, (c) dipoles and octupoles, and
(d) triakontadipoles as a function of $T$, where all multipoles $O_{\Gamma}$
are normalized as ${\rm Tr}[O_{\Gamma}^{2}]=6$ so as to compare all
susceptibilities with the same footing. The susceptibilities only for the TR-
even monopole $I$ and hexadecapole $H_{0}$ [Figs. 2(a) and (b)], the
$z(x)$-components of TR-odd multipoles $J_{z}(J_{x})$,
$T_{z}^{\alpha}(T_{x}^{\alpha})$, $D_{z}^{1\alpha}(D_{x}^{1\alpha})$,
$D_{z}^{2\alpha}(D_{x}^{2\alpha})$ [Figs. 2(c) and (d)] show the Curie-like
behavior in proportion to $1/T$ at low temperatures, since they have finite
matrix elements between the lowest KD. Other susceptibilities such as
$O_{yz/zx},O_{xy}$ and $H_{z}^{\beta},H_{x}^{\beta}$ [Figs. 2(a) and (b)]
increase slightly with decreasing $T$ and then saturate to constant values,
since they only have finite matrix elements between the lowest KD and the
excited KDs, corresponding to Van-Vleck term. As a consequence, active
multipoles at low temperatures are limited to those belonging to the same IRR
under $C_{4v}$: TR-even $\Gamma_{1}$ multipoles $(I,~{}H_{0},~{}O_{u},H_{4})$,
and TR-odd $\Gamma_{2}$ and $\Gamma_{5}$ multipoles,
$(J_{z},~{}T_{z}^{\alpha},~{}D_{z}^{1\alpha,2\alpha})$ and
$(J_{x/y},~{}T_{x/y}^{\alpha},~{}D_{x/y}^{1\alpha,2\alpha})$, respectively.

Figure 2: (Color online) $T$-dependence of $\chi_{O_{{}_{\Gamma}}}^{(0)}$ for
(a) monopole and quadrupoles, (b) hexadecapoles, (c) dipoles and octupoles,
and (d) triakontadipoles, where the $y$-component of
$\chi_{O_{{}_{\Gamma}}}^{(0)}$ (not shown) is the same as the $x$-component
due to the tetragonal symmetry, for example
$\chi_{J_{y}}^{(0)}=\chi_{J_{x}}^{(0)}$.  Figure 3: (Color online)
$\bm{q}$-dependence of the uniform interaction $K_{O_{{}_{\Gamma}}}^{({\rm
u})}(\bm{q})$ for (a),(b) TR-even and (c)-(e) TR-odd multipoles and the
staggered interaction $K_{O_{{}_{\Gamma}}}^{({\rm s})}(\bm{q})$ for (f),(g)
TR-even and (h)-(j) TR-odd multipoles along the high symmetry line in the BZ.

Then we calculate the RKKY interactions between 12 $f$ states/uc (6 CEF states
$m$ and 2 sites $\alpha={\rm Ce1},~{}{\rm Ce2}$) by the fourth-order
perturbation w. r. t. $\mathscr{H}_{c-f}$[11, 12] which is given by,

$\displaystyle\\!\\!K_{m_{1}m_{2}m_{3}m_{4}}^{\alpha,\beta}(\bm{q})\\!=\\!\frac{4}{\Delta_{0}^{2}}\frac{1}{N}\\!\sum_{\bm{k}nn^{\prime}}\\!\mathcal{M}^{n\bm{k}}_{\beta
m_{3},\alpha m_{1}}\mathcal{M}^{n\bm{k}+\bm{q}}_{\alpha m_{2},\beta
m_{4}}\frac{f_{n^{\prime}\bm{k}+\bm{q}}-f_{n\bm{k}}}{\varepsilon_{n\bm{k}}-\varepsilon_{n^{\prime}\bm{k}+\bm{q}}},$
(3)

where $\mathcal{M}^{n\bm{k}}_{\alpha m,\beta
m^{\prime}}=\sum_{\ell\ell^{\prime}}V_{\bm{k}\ell\alpha
m}^{*}V_{\bm{k}\ell^{\prime}\beta
m^{\prime}}U_{n\bm{k}\ell}^{c*}U_{n\bm{k}\ell^{\prime}}^{c}$ and
$V_{\bm{k}\ell\alpha m}$ is the one-body matrix element of
$\mathscr{H}_{c\textrm{-}f}$ between $c~{}(\ell)$ and $f~{}(\alpha m)$ states
and $U_{n\bm{k}\ell}^{c}$ is the eigenvector of $\mathscr{H}_{c\textrm{-}c}$
for $\ell$. $f_{n\bm{k}}$ is the Fermi distribution function
$f_{n\bm{k}}=[\exp\\{(\varepsilon_{n\bm{k}}-\mu)/\tilde{T}\\}+1]^{-1}$, where
we fix the temperature used for Eq. (3) as $\tilde{T}=10~{}{\rm meV}$, since
the RKKY interaction seems less temperature dependent in the present nesting-
free FSs. By using Eq. (3), we obtain the interactions for the uniform
(staggered) multipole $O_{\Gamma}^{({\rm u})}(O_{\Gamma}^{({\rm
s})})=[O_{\Gamma}^{\rm Ce1}+(-)O_{\Gamma}^{\rm Ce2}]/\sqrt{2}$,

$\displaystyle K_{O_{{}_{\Gamma}}}^{({\rm
u,s})}(\bm{q})=\frac{1}{2}\sum_{\alpha\beta}\left[\left\\{\delta_{\alpha\beta}(1\mp
1)\pm
1\right\\}K_{O_{{}_{\Gamma}}}^{\alpha,\beta}(\bm{q})-\delta_{\alpha\beta}K_{O_{{}_{\Gamma}}}^{\alpha,\rm
loc}\right],$ (4)

where the double sign of up/down corresponds to u/s and
$K_{O_{{}_{\Gamma}}}^{\alpha,\beta}(\bm{q})=(1/6)\sum_{\\{m_{i}\\}}(O^{\alpha}_{\Gamma})_{m_{1}m_{2}}K_{m_{1}m_{2}m_{3}m_{4}}^{\alpha,\beta}(\bm{q})(O_{\Gamma}^{\beta\dagger})_{m_{3}m_{4}}$
and $K_{O_{\Gamma}}^{\alpha,\rm
loc}=(1/N)\sum_{\bm{q}}K_{O_{{}_{\Gamma}}}^{\alpha,\alpha}(\bm{q})$. In
calculations, we use $\bm{k}$-mesh of $N=64^{3}$ and consider only $f^{0}$
intermediate excited state, since the contribution of $f^{2}$-process is found
to be not primary for determining the multipole order from the study on CeB6
based on the dynamical mean field theory with Hubbard I approximation[21]. The
excitation energy is set to $\Delta_{0}=2$ eV corresponding to the center of
mass of the $f$ bands of LaCoSi shown in Fig. 1(a), which results in an energy
denominator factor of 1 for the interaction $K(\bm{q})$ (see §4 of SM[16]).

The obtained $K_{O_{{}_{\Gamma}}}^{(\rm u)}(\bm{q})$ and
$K_{O_{{}_{\Gamma}}}^{(\rm s)}(\bm{q})$ for all multipoles listed in Table 1
are plotted along the high symmetry line in BZ shown in Figs. 3(a)-(e) and
(f)-(j), respectively, where the positive (negative) value for a certain
multipole $O_{\Gamma}^{({\rm u,s})}$ with $\bm{q}$ enhances (suppresses) the
corresponding multipole ordering tendency. All $K_{O_{{}_{\Gamma}}}^{({\rm
u,s})}(\bm{q})$ take the values in the range of $-2$ meV to $3$ meV and have
nontrivial $\bm{q}$ and multipole dependence reflecting the $c$ band
$\varepsilon_{n\bm{k}}$ and $c$-$f$ mixing $V_{\bm{k}\ell\alpha m}$ in the
present DFT-based TB model with the symmetry of the tetragonal crystal, for
example $O_{yz}$ and $O_{zx}$ have the same (different) values in the
$\Gamma$-M ($\Gamma$-X) direction. Many multipoles show a maximum at
$\bm{q}=\bm{0}$: the largest component is the staggered dipole $J_{z}$ [Fig.
3(h)], the second largest is the uniform hexadecapole $H_{0}$ [Fig. 3(a)] and
staggered monopole $I$ [Fig. 3(f)], the third is the staggered triakontadipole
$D_{x,y}^{1\alpha}$ [Fig. 3(i)]. These values consist mainly of four nearest-
neighbor couplings between Ce1-Ce2 and four next-nearest-neighbor couplings
between Ce1-Ce1(Ce2-Ce2), whose signs and magnitudes are complicated by the
type of multipoles. The detailed structure of these interactions in real space
and their relation to FSs and $c$-$f$ mixing will be reported somewhere in the
near future.

The RPA susceptibility matrix $\hat{\chi}(\bm{q})$ is given with the matrices
of the single-site susceptibility $\hat{\chi}^{(0)}$ and the interaction
$\hat{K}(\bm{q})$ as follows,

$\displaystyle\hat{\chi}(\bm{q})=\hat{\chi}^{(0)}\left[\hat{1}-\hat{K}(\bm{q})\hat{\chi}^{(0)}\right]^{-1},$
(5)

where each matrix of Eq. (5) has $2\times 6^{2}=72$ dimension corresponding to
$i=(\alpha m_{1}m_{2})$ of the onsite $f$ basis, leading to
$[\hat{\chi}^{(0)}]_{i,j}=\delta_{\alpha\beta}\chi_{m_{1}m_{2}m_{3}m_{4}}^{(0)}$
and
$[\hat{\chi}(\hat{K})]_{i,j}=\chi_{m_{1}m_{2}m_{3}m_{4}}^{\alpha,\beta}~{}(K_{m_{1}m_{2}m_{3}m_{4}}^{\alpha,\beta})$.
By diagonalizing the matrix $\hat{K}(\bm{q})\hat{\chi}^{(0)}$ at each
$\bm{q}$, we obtain the maximum (first largest) eigenvalue
$\alpha_{1}(\bm{q})$ which is enhanced towards a certain ordering instability
with the ordered wavevector $\bm{q}=\bm{q}_{\rm max}$ and finally reaches
$\alpha_{1}=1$ for the critical point of the multipole ordering phase
transition.

Figure 4(a) shows the $i$-th largest eigenvalue $\alpha_{i}(\bm{q}_{\rm max})$
up to 3rd values as a function of $T$. Each eigenvalue increases with
decreasing $T$ and the maximum eigenvalue $\alpha_{1}=1$ is realized at
$T=0.7$ meV with $\bm{q}_{\rm max}=\bm{0}$. Here, $T$-dependence of
$\alpha_{i}(\bm{q})$ is due to the single-site susceptibility $\chi^{(0)}$
shown in Fig. 2, since the interaction $K(\bm{q})$ is independent of
temperature as mentioned above. The interactions of the active multipoles at
low temperatures show maximums at $\bm{q}=\bm{0}$ as shown in Fig. 3, so that
the corresponding multipole fluctuations develop towards the multipole order
at $\bm{q}=\bm{0}$. For explicit check, we plot the $\bm{q}$-dependence of
$\alpha_{i}(\bm{q})$ shown in the inset of Fig. 4(b) for
$(\alpha_{1},~{}T)=(0.95,~{}0.73~{}{\rm meV})$, where the clear peaks can be
seen at $\bm{q}=\bm{0}$ for each $\alpha_{i}$ and the first and second
eigenvalues are quite close to each other.

Figure 4: (Color online) (a) $T$-dependence of the $i$-th largest eigenvalue
$\alpha_{i}(\bm{q}_{\rm max})$ of the matrix $\hat{K}(\bm{q})\hat{\chi}^{(0)}$
with $\bm{q}_{\rm max}=\bm{0}$. (c),(e) $T$-dependence of the enhancement
rates of the $\bm{q}=\bm{0}$ RPA susceptibility for (c) uniform and (e)
staggered multipoles $(I,~{}H_{0},~{}O_{u},~{}H_{4})$ of the TR-even
$\Gamma_{1}$. (b),(d),(f) The corresponding $\bm{q}$-dependences for
$(\alpha_{1},T)=(0.95,~{}0.73~{}{\rm meV})$.

In order to clarify such $\bm{q}=\bm{0}$ order obtained in Fig. 4(a), we
examine the enhancement rate of the RPA susceptibility
$\chi_{O_{{}_{\Gamma}}}^{({\rm u,s})}(\bm{q})/\chi_{O_{{}_{\Gamma}}}^{(0)}$
for all multipoles. As a results, we have found that the RPA susceptibilities
for TR-even $\Gamma_{1}$ multipoles $(I,~{}H_{0},~{}O_{u},~{}H_{4})$ dominate
over the other multipoles.

Figures 4(c)-(f) show the enhancement rates for the TR-even $\Gamma_{1}$
multipoles of (c) uniform and (e) staggered components as a function of $T$,
where Figs. 4(d) and (f) show the $\bm{q}$-dependence of the corresponding
enhancement rates for the same values of Fig. 4(b). With decreasing $T$,
$\chi_{O_{{}_{\Gamma}}}^{({\rm u})}/\chi_{O_{{}_{\Gamma}}}^{(0)}$ of the
monopole $I$ and hexadecapole $H_{0}$ increase rapidly with keeping the same
value as shown in Fig. 4(c). A similar enhancement is also observed in
$\chi_{O_{{}_{\Gamma}}}^{({\rm s})}/\chi_{O_{{}_{\Gamma}}}^{(0)}$ for $I$ and
$H_{0}$ as shown in Fig. 4(e), but the staggered component is larger than the
uniform one as compered to Fig. 4(c) and (d), where $\chi_{I}^{({\rm
s})}/\chi_{I}^{(0)}\simeq\chi_{H_{0}}^{({\rm s})}/\chi_{H_{0}}^{(0)}\simeq 20$
shown in Fig. 4(f) is consistent with the enhancement rate from Eq. (5)
$1/(1-\alpha_{1})=20$ times, corresponding to the maximum eigenmode in Fig.
4(b). As a consequence, an antiferro multipole order mixed with $(I,~{}H_{0})$
is realized due to the $\bm{q}=\bm{0}$ staggered fluctuation as the maximum
RPA mode, while the $\bm{q}=\bm{0}$ uniform fluctuation of $(I,~{}H_{0})$
inducing a phase separation instability is the second largest one which is not
realized in actually. The fluctuations of $O_{u}$ and $H_{4}$ also increase
with keeping slightly small value compared to the main components $I$ and
$H_{0}$ as shown in Fig. 4(c)-(f), since the single-site susceptibility and
interactions for different multipoles belonging to the same IRR $\Gamma_{1}$
are also finite (not shown). We also investigate the rates for other
multipoles except for the TR-even $\Gamma_{1}$ explicitly shown in §5 of
SM[16], where the rates for the quadrupoles and hexadecapoles have a small
$T$-dependence and are hardly enhanced, while that for TR-odd multipoles of
$(J_{z},~{}T_{z}^{\alpha},~{}D_{z}^{1\alpha,2\alpha})~{}[\Gamma_{2}]$ and
$(J_{x/y},~{}T_{x/y}^{\alpha},~{}D_{x/y}^{1\alpha,2\alpha})~{}[\Gamma_{5}]$
enhance two to three times.

Figure 5: (Color online) Schematic pictures of the $f$-CDW order in the (a)
tetragonal and (b) triclinic lattices, where $f$-electron electric dipoles
$\bm{p}_{f}$ are dipicted by (a) dotted and (b) solid vectors, and large red
(small blue) circles indicate an increase (decrease) in the number of $f$
electrons at Ce atoms.

Here we discuss the relationship between the present results and previous
studies and experiments. First, the obtained antiferro order of $(I,~{}H_{0})$
is consistent with the HO in the experiments in the sense of nonmagnetic[1, 2,
3]. This $(I,~{}H_{0})$ order is also a $f$-electron charge density wave
($f$-CDW) with $\bm{q}=\bm{0}$, yielding the $f$-electron number $n_{f}$
imbalance on the two Ce atoms in uc, where $n_{f}=1$ on each Ce atom in
$T>T_{0}$ splits into $n_{f}=1+(-)\delta$ on the Ce1 (Ce2) atom in $T<T_{0}$
as shown in Fig. 5(a). Similar antiferro charge ordering in $f$-electron
systems has been discussed in studies of the filled skutterudite compounds
PrRu4P12[22] and SmRu4P12[23], where $f$-electron charge ordering is driven by
the perfect nesting of $c$-electron FS and realized in both metal-insulator
transition. This is contrast to the present $f$-CDW driven by the RKKY
interactions derived from the DFT bandstructure without the nesting, where
$\bm{q}=\bm{0}$ charge ordering is realized on the two Ce atoms in uc and the
$c$-electron system remains metal.

The deviation of $n_{f}$ from 1 in the present $f$-CDW is considered to be
small ($\delta\ll 1$), and some of them are expected to transfer to the $c$
electron system, i. e., the $f$ electrons should become somewhat itinerant for
$T<T_{0}$, which is supported by the optical conductivity experiment[24]. In
addition, the present $f$-CDW is coupled to the hexadecapole $H_{0}$, yielding
to a shift in the CEF levels with an opposite splitting at the two Ce atoms,
which may be related to the mysterious excitation observed in neutron
scattering[7]. Moreover, the contributions of $I$ and $H_{0}$ to the $f$-CDW
would vary gradually with pressure, which could also explain the $T_{0}$ dip
observed in electrical resistivity under high pressure[25] by accompanying the
CDW of the $c$-electrons, where the $f$-electrons are more itinerant.

Finally, under the present $f$-CDW order, four equivalent electric dipoles
$\bm{p}_{f}$ connecting Ce1-Ce2 are realized as shown in Fig. 5(a), and when
the triclinic distortion occurs one of them become inequivalent and a single
electric dipole with a vector connecting the intracell Ce atoms in uc is
chosen as shown in Fig. 5(b). The triclinic structural phase transition in the
experiments[8, 9] may be regarded as a consequence of the energy stabilizing
of such the induced electric dipole alignment. Alternatively, it may be
necessary to introduce the electron-lattice interactions into the present
model in order to describe the structural transition and is one of the
important issues to be addressed in the future. A unified understanding of
such the structural phase transition at ambient pressure and anomalies at high
pressure is an important issue that can be positioned as an extension of this
study.

In summary, we have studied the RKKY interactions between $f$ electrons based
on the realistic Wannier model derived from the ab initio calculations on
CeCoSi. The RPA analysis reveals that the leading order is a staggered
$\bm{q}=\bm{0}$ monopole order coupled with hexadecapole $H_{0}$, which
corresponds to an $f$-CDW with a staggered shift in the crystal field
splitting and can explain the some experiments in the HO.

We would like to thank H. Tanida and for valuable comments and fruitful
discussions. This work was supported by JSPS KAKENHI Grant Numbers JP24K06943
and JP21H01031.

## References

  * [1] H. Tanida, Y. Muro, and T. Matsumura, J. Phys. Soc. Jpn. 87, 023705 (2018).
  * [2] H. Tanida, K. Mitsumoto, Y. Muro, T. Fukuhara, Y. Kawamura, A. Kondo, K. Kindo, Y. Matsumoto, T. Namiki, T. Kuwai, and T. Matsumura, J. Phys. Soc. Jpn. 88, 054716 (2019).
  * [3] M. Manago, H. Kotegawa, H. Tou, H. Harima, and H. Tanida, J. Phys. Soc. Jpn. 90, 023702 (2021).
  * [4] M. Yatsushiro and S. Hayami, J. Phys. Soc. Jpn. 89, 013703 (2020).
  * [5] M. Yatsushiro and S. Hayami, Phys. Rev. B 102, 195147 (2020).
  * [6] H. Hidaka, S. Yanagiya, E. Hayasaka, Y. Kaneko, T. Yanagisawa, H. Tanida, and H. Amitsuka, J. Phys. Soc. Jpn. 91, 094701 (2022).
  * [7] S. E. Nikitin, D. G. Franco, J. Kwon, R. Bewley, A. Podlesnyak, A. Hoser, M. M. Koza, C. Geibel, and O. Stockert, Phys. Rev. B 101, 214426 (2020).
  * [8] T. Matsumura, S. Kishida, M. Tsukagoshi, Y. Kawamura, H. Nakao, and H. Tanida, J. Phys. Soc. Jpn. 91, 064704 (2022).
  * [9] M. Manago, A. Ishigaki, H. Tou, H. Harima, H. Tanida, and H. Kotegawa, Phys. Rev. B 108, 085118 (2023).
  * [10] M. Yatsushiro and S. Hayami, J. Phys. Soc. Jpn. 91, 104701 (2022).
  * [11] T. Yamada and K. Hanzawa, J. Phys. Soc. Jpn. 88, 084703 (2019).
  * [12] T. Yamada and K. Hanzawa, JPC Conf. 30, 011150 (2020).
  * [13] M. Kawamura, Comput. Phys. Commun. 239, 197 (2019).
  * [14] P. Blaha, K. Schwarz, F. Tran, R. Laskowski, G. K. H. Madsen, and L. D. Marks, J. Chem. Phys. 152, 074101 (2020).
  * [15] G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. Mostofi, and J. Yates, J. Phys.: Condens. Mater 32, 165902 (2020).
  * [16] (Supplemental Material) The technical details of the calculations and supplemental results are available online(see, also, Refs. [[26, 27]] therein).
  * [17] A. Kimura, T. Yoshikawa, Y. Fukushima, K. Shiraishi, T. Kono, M. Nurmamat, M. Arita, T. Muro, K.Mitsumoto, and H. Tanida, SCES2020 -International Conference on Strongly Correlated Electron Systems, oral, F4C (2021).
  * [18] H. Kusunose, J. Phys. Soc. Jpn. 77, 064710 (2008).
  * [19] S. Hayami, M. Yatsushiro, Y. Yanagi, and H. Kusunose, Phys. Rev. B 98, 165110 (2018).
  * [20] H. Ikeda, M.-T. Suzuki, R. Arita, T. Takimoto, T. Shibauchi, and Y. Matsuda, Nat. Phys. 8, 528 (2012).
  * [21] J. Otsuki, K. Yoshimi, H. Shinaoka, and O. Jeschke, arXiv:2209.10429 (2022).
  * [22] T. Takimoto, J. Phys. Soc. Jpn. 75, 034714 (2006).
  * [23] R. Shiina, J. Phys. Soc. Jpn. 82, 083713 (2013).
  * [24] S. Kimura, H. Watanabe, S. Tatsukawa, and H. Tanida, J. Phys. Soc. Jpn. 92, 043704 (2023).
  * [25] E. Lengyel, M. Nicklas, N. Caroca-Canales, C. G. F. Weng, M. Smidman, L. Jiao, X. Lu, and H. Q. Yuan, Phys. Rev. B 88, 155137 (2013).
  * [26] J. P. Perdew, S. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * [27] K. Hanzawa, J. Phys. Soc. Jpn. 84, 024717 (2015).

Supplemental material to “RKKY Interactions and Multipole Order in Ab initio
Wannier Model of CeCoSi”
Takemi Yamada, Yuki Yanagi and Keisuke Mitsumoto

Liberal Arts and Sciences, Toyama Prefectural University, Imizu, Toyama
939-0398, Japan

This supplemental material contains: 1\. DFT calculation conditions for
WIEN2k, 2\. TB model construction for Wannier90, 3\. Multipole operator and
single-site susceptibility, 4\. Derivation of RKKY interactions and 5\.
Enhancement rates of the RPA susceptibilities.

##  1 DFT calculation for WIEN2k

Here we provide the details of the DFT calculation based on WIEN2k [14], which
is the all-electron first-principles code, where the basis functions are
expanded by the relativistic full-potential augmented plane wave (FLAPW)
and/or APW + local orbitals depending on the core and valence states of each
atom in uc. We employ the generalized gradient approximation for the exchange-
correlation potential of the PBE-GGA potential [26]. The SOC is fully included
within the second variation approximation.

In self-consistent calculation, we use the experimental lattice parameters[2]
in the space group $P4/nmm$ (#129), the lattice constants
$a=b=4.057~{}{\rm\AA}$ and $c=6.987~{}{\rm\AA}$, the internal fractional
coordinates $(x,y,z)=(0.25,0.25,0.6781)$ for Ce, $(0.75,0.25,0)$ for Co and
$(0.25,0.25,0.183)$ for Si, and the calculation conditions: $\bm{k}$-mesh of
$20^{2}\times 23$ corresponding to 2772 $\bm{k}$-points in the irreducible
part of BZ and the muffin-tin radius $R_{\rm MT}^{\rm Ce,Co(Si)}=2.5~{}(1.91)$
bohr. The plane wave cutoff is $R_{\rm MT}K_{\rm max}=9$ and the energy and
charge convergence conditions are set to 10-6 Ry and 10-5 e, respectively.

##  2 TB model construction by Wannier90

We use Wannier90[15] to construct a TB model of CeCoSi with localized Wannier
orbitals, which are given by the inverse Fourier transform of the DFT Bloch
states and thus can fully reproduce the DFT bands within a certain energy
window. In particular, well localized Wannier orbitals at each atomic center
can be treated as atomic orbitals, allowing for more detailed microscopic
analysis than the DFT band calculation itself and incorporating many-body
interactions not included in the DFT calculation, such as the RKKY interaction
calculated in this study.

In order to simultaneously achieve the main valence bandstructure and Wannier
function localization at a high level, we construct a TB model with 84
localized Wannier orbitals, which consists of Ce-$4f,5d~{}(7+5=12)$,
Co-$3d~{}(5)$, and Si-$3p,4s~{}(3+1=4)$ for a total of $12+5+4=21$ orbitals
with spin (2) degrees of freedom and the sub-lattice (2) degrees of freedom by
2CeCoSi/uc, resulting in $21\times 2\times 2=84$ orbitals. The inner window
(dis-frozen energy window) corresponding to the region reproducing the
bandstructure is set in the energy range from $-10~{}{\rm eV}$ to $+4~{}{\rm
eV}$ and the target band for adjusting the outer window is set to
$84+180=264$, where the localization of the Wannier function improves as the
target band increases, though in this study the qualitative results are almost
unchanged by the case of about $84+20=104$ bands. The $\bm{k}$-mesh is set to
$8^{2}\times 6=384$ corresponding to $567$ real-space hopping vectors and the
disentangle convergence condition is set to 10${}^{-7}~{}\AA^{2}$.

The obtained TB model is given by the following form,

$\displaystyle\mathscr{H}_{\rm
TB}=\mathscr{H}_{c\textrm{-}c}+\mathscr{H}_{f\textrm{-}f}+\mathscr{H}_{c\textrm{-}f},$
(6)
$\displaystyle\quad\mathscr{H}_{c\textrm{-}c}=\sum_{\bm{i},\bm{\delta}}\sum_{\ell\ell^{\prime}}h_{\bm{i}\ell\ell^{\prime}}^{cc}(\bm{\delta})c_{\bm{i}\ell}^{\dagger}c_{\bm{i}+\bm{\delta}\ell^{\prime}},$
(7)
$\displaystyle\quad\mathscr{H}_{f\textrm{-}f}=\sum_{\bm{i},\bm{\delta}}\sum_{mm^{\prime}}h_{\bm{i}mm^{\prime}}^{ff}(\bm{\delta})f_{\bm{i}m}^{\dagger}f_{\bm{i}+\bm{\delta}m^{\prime}},$
(8)
$\displaystyle\quad\mathscr{H}_{c\textrm{-}f}=\sum_{\bm{i},\bm{\delta}}\sum_{\ell
m}\left(V_{\bm{i}\ell
m}(\bm{\delta})c_{\bm{i}\ell}^{\dagger}f_{\bm{i}+\bm{\delta}m}+h.c.\right),$
(9)

where $c_{\bm{i}\ell}^{\dagger}~{}(f_{\bm{i}m}^{\dagger})$ is a creation
operator of $c~{}(f)$ electron in uc at $\bm{R_{i}}$ with a state $\ell~{}(m)$
being 56 $d,p,s$ real functions per spin on each atom (28 CEF eigenstates with
total angular momentum $j=5/2,7/2$ on each Ce atom). The $c$-$c$ ($f$-$f$)
matrix element of
$h_{\bm{i}\ell\ell^{\prime}}^{cc}~{}(h_{\bm{i}mm^{\prime}}^{ff})$ includes the
$c~{}(f)$ energy levels, SOC couplings, CEF splittings and $c$-$c$ ($f$-$f$)
hopping integrals with the real-space hopping vector $\bm{\delta}$.
$V_{\bm{i}\ell m}(\bm{\delta})$ is the $c$-$f$ mixing element, which is finite
in both the onsite and intersite components due to the lack of inversion
symmetry at each atom.

The wavevector $\bm{k}$-representation of $\mathscr{H}_{\rm TB}$ is given by,

$\displaystyle\mathscr{H}_{\rm TB}$
$\displaystyle=\sum_{\bm{k}}\sum_{\ell\ell^{\prime}}h_{\ell\ell^{\prime}}^{cc}(\bm{k})c_{\bm{k}\ell}^{\dagger}c_{\bm{k}\ell^{\prime}}+\sum_{\bm{k}}\sum_{mm^{\prime}}h_{mm^{\prime}}^{ff}(\bm{k})f_{\bm{k}m}^{\dagger}f_{\bm{k}m^{\prime}}+\sum_{\bm{k}}\sum_{\ell
m}\left(V_{\bm{k}\ell m}c_{\bm{k}\ell}^{\dagger}f_{\bm{k}m}+h.c.\right),$ (10)

and $\mathscr{H}_{\rm TB}$ can be diagonalized in $\bm{k}$-space and
transformed into the eigenstates with band-index $n$ as follows,

$\displaystyle\mathscr{H}_{\rm TB}$
$\displaystyle=\sum_{n\bm{k}}E_{n\bm{k}}a_{n\bm{k}}^{\dagger}a_{n\bm{k}},\quad
a_{n\bm{k}}=\sum_{m}U_{n\bm{k}m}f_{\bm{k}m}+\sum_{\ell}U_{n\bm{k}\ell}c_{\bm{k}\ell},$
(11)

where the eigenenergy $E_{n\bm{k}}$ is the full TB bands reproducing the DFT
bands of CeCoSi and an annihilation operator $a_{n\bm{k}}$ can be written by
the the eigenvector components $U_{n\bm{k}m}$ and $U_{n\bm{k}\ell}$. Here we
note that $U_{n\bm{k}\ell}$ is the eigenvector of $E_{n\bm{k}}$ and different
from $U_{n\bm{k}\ell}^{c}$ which is the eigenvector of $\varepsilon_{n\bm{k}}$
in the main text.

The obtained $E_{n\bm{k}}$ describe the itinerant $f$-electron bands and
disagree with the electronic structure suggested by the ARPES experiment[17],
while as mentioned in the main text, the diagonalized $c$-electron bands of
only $\mathscr{H}_{c\textrm{-}c}$ shows relatively good correspondence with
the ARPES[17], which justifies starting from the localized $f$-electron limit
where $f$ electrons are well localized at each Ce.

##  3 Multipole operator and single-site susceptibility

Here we briefly address the multipole operator and single-site susceptibility
described in the main text. As multipoles in this study, we treat only the
traditional local multipoles given by $6\times 6$ matrices for a Ce-$4f$ one-
electron basis $|j=5/2,m=-j,\cdots,j-1,j\rangle$, i.e., electric multipoles
and magnetic multipoles[18, 19]. In this case, we can define $6\times 6=36$
independent operators described by a spherical tensor operator $T_{q}^{(k)}$
with the same symmetry as the spherical harmonic function on isolated ions in
a vacuum, where $k$ is rank and $q=-k,-k+1,\cdots,+k$. In crystals, we can
also define the tesseral tensor operator $T_{kq}^{\rm(c,s)}$ given by a linear
combination of $T_{q}^{(k)}$ as follows,

$\displaystyle T_{k0}^{\rm(c)}$ $\displaystyle=T_{0}^{(k)},$ (12)
$\displaystyle T_{kq}^{\rm(c)}$
$\displaystyle=(-1)^{q}\frac{1}{\sqrt{2}}\left\\{\left(T_{q}^{(k)}\right)^{\dagger}+T_{q}^{(k)}\right\\}=\begin{cases}\cfrac{1}{\sqrt{2}}\left(T_{-q}^{(k)}-T_{q}^{(k)}\right)&(q={\rm
odd})\\\ \cfrac{1}{\sqrt{2}}\left(T_{-q}^{(k)}+T_{q}^{(k)}\right)&(q={\rm
even})\\\ \end{cases},$ (13) $\displaystyle T_{kq}^{\rm(s)}$
$\displaystyle=(-1)^{q}\frac{i}{\sqrt{2}}\left\\{\left(T_{q}^{(k)}\right)^{\dagger}-T_{q}^{(k)}\right\\}=\begin{cases}\cfrac{i}{\sqrt{2}}\left(T_{-q}^{(k)}+T_{q}^{(k)}\right)&(q={\rm
odd})\\\ \cfrac{i}{\sqrt{2}}\left(T_{-q}^{(k)}-T_{q}^{(k)}\right)&(q={\rm
even})\\\ \end{cases},$ (14)

where $\left(T_{q}^{(k)}\right)^{\dagger}=T_{-q}^{(k)}$ and the linear
combination of the even (odd) rank multipoles corresponds to the electric
(magnetic) multipoles.

Table 2 summarizes the classification of multipole operators in CeCoSi under
the IRR of the point group $C_{4v}$ at each Ce site and their explicit
representation by tesseral tensor operators, where the multipole operators are
written in the same notation as the study on the similar tetragonal compound
URu2Si2, such as monopole (dipole), quadrupole (octupole) and hexadecapole
(triakontadipole) are denoted by $I~{}(J)$, $O~{}(T)$, and $H~{}(D)$,
respectively[20].

Table 2:  List of all multipole operators $O_{\Gamma}$ up to rank 5 and their
IRRs in the point groups $O$ and $C_{4v}$ with $(x,y,z)$\- and
$T_{kq}^{\rm(c,s)}$-notations, where the double sign of IRR
$\Gamma_{\gamma}^{\pm}$ corresponds to the TR-even/TR-odd symmetry.

rank | $O$ | $C_{4v}$ | $O_{\Gamma}$ |  $(x,y,z)$-notation |  $T_{kq}^{\rm(c,s)}$-notation
---|---|---|---|---|---
0 | $\Gamma_{1}^{+}(A_{1}^{+})$ | $\Gamma_{1}^{+}(A^{+})$ | $I$ | $1$ |  $T_{00}^{\rm(c)}$
1 | $\Gamma_{4}^{-}(T_{1}^{-})$ | $\Gamma_{5}^{-}(E^{-})$ | $J_{x}$ | $x$ |  $T_{11}^{\rm(c)}$
|  | $\Gamma_{5}^{-}(E^{-})$ | $J_{y}$ | $y$ |  $T_{11}^{\rm(s)}$
|  | $\Gamma_{2}^{-}(A_{2}^{-})$ | $J_{z}$ | $z$ |  $T_{10}^{\rm(c)}$
2 | $\Gamma_{3}^{+}(E^{+})$ | $\Gamma_{1}^{+}(A_{1}^{+})$ | $O_{u}$ | $\tfrac{1}{2}(3z^{2}-r^{2})$ |  $T_{20}^{\rm(c)}$
|  | $\Gamma_{3}^{+}(B_{1}^{+})$ | $O_{v}$ | $\tfrac{\sqrt{3}}{2}(x^{2}-y^{2})$ |  $T_{22}^{\rm(c)}$
| $\Gamma_{5}^{+}(T_{2}^{+})$ | $\Gamma_{5}^{+}(E^{+})$ | $O_{yz}$ | $\sqrt{3}yz$ |  $T_{21}^{\rm(s)}$
|  | $\Gamma_{5}^{+}(E^{+})$ | $O_{zx}$ | $\sqrt{3}zx$ |  $T_{21}^{\rm(c)}$
|  | $\Gamma_{4}^{+}(B_{2}^{+})$ | $O_{xy}$ | $\sqrt{3}xy$ |  $T_{22}^{\rm(s)}$
3 | $\Gamma_{2}^{-}(A_{2}^{-})$ | $\Gamma_{3}^{-}(A_{2}^{-})$ | $T_{xyz}$ | $\sqrt{15}xyz$ |  $T_{32}^{\rm(s)}$
| $\Gamma_{4}^{-}(T_{1}^{-})$ | $\Gamma_{5}^{-}(E^{-})$ | $T_{x}^{\alpha}$ | $\tfrac{1}{2}x(5x^{2}-3r^{2})$ |  $\sqrt{\tfrac{5}{8}}T_{33}^{\rm(c)}-\sqrt{\tfrac{3}{8}}T_{31}^{\rm(c)}$
|  | $\Gamma_{5}^{-}(E^{-})$ | $T_{y}^{\alpha}$ | $\tfrac{1}{2}y(5y^{2}-3r^{2})$ |  $-\sqrt{\tfrac{5}{8}}T_{33}^{\rm(s)}-\sqrt{\tfrac{3}{8}}T_{31}^{\rm(s)}$
|  | $\Gamma_{2}^{-}(A_{2}^{-})$ | $T_{z}^{\alpha}$ | $\tfrac{1}{2}z(5z^{2}-3r^{2})$ |  $T_{30}$
| $\Gamma_{5}^{-}(T_{2}^{-})$ | $\Gamma_{5}^{-}(E^{-})$ | $T_{x}^{\beta}$ | $\tfrac{\sqrt{15}}{2}x(y^{2}-z^{2})$ |  $-\sqrt{\tfrac{3}{8}}T_{33}^{\rm(c)}-\sqrt{\tfrac{5}{8}}T_{31}^{\rm(c)}$
|  | $\Gamma_{5}^{-}(E^{-})$ | $T_{y}^{\beta}$ | $\tfrac{\sqrt{15}}{2}y(z^{2}-x^{2})$ |  $-\sqrt{\tfrac{3}{8}}T_{33}^{\rm(s)}+\sqrt{\tfrac{5}{8}}T_{31}^{\rm(s)}$
|  | $\Gamma_{4}^{-}(B_{2}^{-})$ | $T_{z}^{\beta}$ | $\tfrac{\sqrt{15}}{2}z(x^{2}-y^{2})$ |  $T_{32}^{\rm(c)}$
4 | $\Gamma_{1}^{+}(A_{1}^{+})$ | $\Gamma_{1}^{+}(A_{1}^{+})$ | $H_{0}$ | $\tfrac{5\sqrt{21}}{12}(x^{4}+y^{4}+z^{4}-\tfrac{3}{5}r^{4})$ |  $\sqrt{\tfrac{5}{12}}T_{44}^{\rm(c)}+\sqrt{\tfrac{7}{12}}T_{40}^{\rm(c)}$
| $\Gamma_{3}^{+}(E^{+})$ | $\Gamma_{1}^{+}(A_{1}^{+})$ | $H_{4}$ | $\tfrac{7\sqrt{15}}{12}[2z^{4}-x^{4}-y^{4}-\tfrac{6}{7}r^{2}(3z^{2}-r^{2})]$ |  $-\sqrt{\tfrac{7}{12}}T_{44}^{\rm(c)}+\sqrt{\tfrac{5}{12}}T_{40}^{\rm(c)}$
|  | $\Gamma_{3}^{+}(B_{1}^{+})$ | $H_{2}$ | $\tfrac{7\sqrt{5}}{4}[x^{4}-y^{4}-\tfrac{6}{7}r^{2}(x^{2}-y^{2})]$ |  $-T_{42}^{\rm(c)}$
| $\Gamma_{4}^{+}(T_{1}^{+})$ | $\Gamma_{5}^{+}(E^{+})$ | $H_{x}^{\alpha}$ | $\tfrac{\sqrt{35}}{2}yz(y^{2}-z^{2})$ |  $-\sqrt{\tfrac{1}{8}}T_{43}^{\rm(s)}-\sqrt{\tfrac{7}{8}}T_{41}^{\rm(s)}$
|  | $\Gamma_{5}^{+}(E^{+})$ | $H_{y}^{\alpha}$ | $\tfrac{\sqrt{35}}{2}zx(z^{2}-x^{2})$ |  $-\sqrt{\tfrac{1}{8}}T_{43}^{\rm(c)}+\sqrt{\tfrac{7}{8}}T_{41}^{\rm(c)}$
|  | $\Gamma_{2}^{+}(A_{2}^{+})$ | $H_{z}^{\alpha}$ | $\tfrac{\sqrt{35}}{2}xy(x^{2}-y^{2})$ |  $T_{44}^{\rm(s)}$
| $\Gamma_{5}^{+}(T_{2}^{+})$ | $\Gamma_{5}^{+}(E^{+})$ | $H_{x}^{\beta}$ | $\tfrac{\sqrt{5}}{2}yz(7x^{2}-r^{2})$ |  $\sqrt{\tfrac{7}{8}}T_{43}^{\rm(s)}-\sqrt{\tfrac{1}{8}}T_{41}^{\rm(s)}$
|  | $\Gamma_{5}^{+}(E^{+})$ | $H_{y}^{\beta}$ | $\tfrac{\sqrt{5}}{2}zx(7y^{2}-r^{2})$ |  $-\sqrt{\tfrac{7}{8}}T_{43}^{\rm(c)}-\sqrt{\tfrac{1}{8}}T_{41}^{\rm(c)}$
|  | $\Gamma_{4}^{+}(B_{2}^{+})$ | $H_{z}^{\beta}$ | $\tfrac{\sqrt{5}}{2}xy(7z^{2}-r^{2})$ |  $T_{42}^{\rm(s)}$
5 | $\Gamma_{3}^{-}(E^{-})$ | $\Gamma_{1}^{-}(A_{1}^{-})$ | $D_{4}$ | $\tfrac{3\sqrt{35}}{2}xyz(x^{2}-y^{2})$ |  $T_{54}^{\rm(s)}$
|  | $\Gamma_{3}^{-}(B_{1}^{-})$ | $D_{2}$ | $-\tfrac{\sqrt{105}}{2}xyz(3z^{2}-r^{2})$ |  $-T_{52}^{\rm(s)}$
| $\Gamma_{4}^{-}(T_{1}^{-})$ | $\Gamma_{5}^{-}(E^{-})$ | $D_{x}^{1\alpha}$ | $\tfrac{x}{8}\left[8x^{4}-40x^{2}(y^{2}+z^{2})+15(y^{2}+z^{2})^{2}\right]$ |  $\tfrac{1}{8\sqrt{2}}(3\sqrt{7}T_{55}^{\rm(c)}-\sqrt{35}T_{53}^{\rm(c)}+\sqrt{30}T_{51}^{\rm(c)})$
|  | $\Gamma_{5}^{-}(E^{-})$ | $D_{y}^{1\alpha}$ | $\tfrac{y}{8}\left[8y^{4}-40y^{2}(z^{2}+x^{2})+15(z^{2}+x^{2})^{2}\right]$ |  $\tfrac{1}{8\sqrt{2}}(3\sqrt{7}T_{55}^{\rm(s)}+\sqrt{35}T_{53}^{\rm(s)}+\sqrt{30}T_{51}^{\rm(s)})$
|  | $\Gamma_{2}^{-}(A_{2}^{-})$ | $D_{z}^{1\alpha}$ | $\tfrac{z}{8}\left[8z^{4}-40z^{2}(x^{2}+y^{2})+15(x^{2}+y^{2})^{2}\right]$ |  $T_{50}$
| $\Gamma_{4}^{-}(T_{1}^{-})$ | $\Gamma_{5}^{-}(E^{-})$ | $D_{x}^{2\alpha}$ | $\tfrac{3\sqrt{35}}{2}x\left[y^{4}+z^{4}-\tfrac{3}{4}(y^{2}+z^{2})^{2}\right]$ |  $\tfrac{1}{16}(\sqrt{10}T_{55}^{\rm(c)}+9\sqrt{2}T_{53}^{\rm(c)}+2\sqrt{21}T_{51}^{\rm(c)})$
|  | $\Gamma_{5}^{-}(E^{-})$ | $D_{y}^{2\alpha}$ | $\tfrac{3\sqrt{35}}{2}y\left[z^{4}+x^{4}-\tfrac{3}{4}(z^{2}+x^{2})^{2}\right]$ |  $\tfrac{1}{16}(\sqrt{10}T_{55}^{\rm(s)}-9\sqrt{2}T_{53}^{\rm(s)}+2\sqrt{21}T_{51}^{\rm(s)})$
|  | $\Gamma_{2}^{-}(A_{2}^{-})$ | $D_{z}^{2\alpha}$ | $\tfrac{3\sqrt{35}}{2}z\left[x^{4}+y^{4}-\tfrac{3}{4}(x^{2}+y^{2})^{2}\right]$ |  $T_{54}^{\rm(c)}$
| $\Gamma_{5}^{-}(T_{2}^{-})$ | $\Gamma_{5}^{-}(E^{-})$ | $D_{x}^{\beta}$ | $\tfrac{\sqrt{105}}{4}x(y^{2}-z^{2})(3x^{2}-r^{2})$ |  $\tfrac{1}{4\sqrt{2}}(-\sqrt{15}T_{55}^{\rm(c)}-\sqrt{3}T_{53}^{\rm(c)}+\sqrt{14}T_{51}^{\rm(c)})$
|  | $\Gamma_{5}^{-}(E^{-})$ | $D_{y}^{\beta}$ | $\tfrac{\sqrt{105}}{4}y(z^{2}-x^{2})(3y^{2}-r^{2})$ |  $\tfrac{1}{4\sqrt{2}}(\sqrt{15}T_{55}^{\rm(s)}-\sqrt{3}T_{53}^{\rm(s)}-\sqrt{14}T_{51}^{\rm(s)})$
|  | $\Gamma_{4}^{-}(B_{2}^{-})$ | $D_{z}^{\beta}$ | $\tfrac{\sqrt{105}}{4}z(x^{2}-y^{2})(3z^{2}-r^{2})$ |  $T_{52}^{\rm(c)}$

Figure 6: (Color online) $T$-dependence of (a) the particle-number
$n_{i}=n_{i+}+n_{i-}$ for three KDs $(i=1,2,3)$ and (b) the functions
$f_{1}$-$f_{6}$ in the single-site susceptibility
$\chi_{m_{1}m_{2}m_{3}m_{4}}^{(0)}$.

The single-site susceptibility $\chi_{m_{1}m_{2}m_{3}m_{4}}^{(0)}$ for 6 CEF
states in the experiment[7] is given by the Lehmann representation as follows,

$\displaystyle\chi_{m_{1}m_{2}m_{3}m_{4}}^{(0)}$
$\displaystyle=\frac{1}{Z}\sum_{mm^{\prime}}A_{mm^{\prime}}^{m_{1}m_{2}}A_{m^{\prime}m}^{m_{4}m_{3}}\left\\{\frac{\delta_{E_{m},E_{m^{\prime}}}}{T}e^{-E_{m}/T}+\left(1-\delta_{E_{m},E_{m^{\prime}}}\right)\frac{e^{-E_{m^{\prime}}/T}-e^{-E_{m}/T}}{E_{m}-E_{m^{\prime}}}\right\\},$
(15)

where $A_{mm^{\prime}}^{m_{1}m_{2}}=\langle
m|f_{m_{1}}^{\dagger}f_{m_{2}}|m^{\prime}\rangle$ and $f_{m}^{\dagger}$ is a
creation operator for a $f$ electron of CEF state $m$ with energy $E_{m}$. $Z$
is a partition function $Z=2(1+e^{-\Delta_{1}/T}+e^{-\Delta_{2}/T})$ with the
1st (2nd) excited level $\Delta_{1(2)}=10.78~{}(14.26)$ meV. The 6 CEF states
of 3 KDs for $\\{m,m^{\prime},m_{i}\\}$ are explicitly given by,

$\displaystyle\Ket{1\pm}=\Ket{\Gamma_{7a\pm}}=\mp\sqrt{1-w^{2}}\ket{\pm\tfrac{5}{2}}\pm
w\ket{\mp\tfrac{3}{2}},$ (16)
$\displaystyle\Ket{2\pm}=\Ket{\Gamma_{7b\pm}}=w\ket{\pm\tfrac{5}{2}}+\sqrt{1-w^{2}}\ket{\mp\tfrac{3}{2}},$
(17) $\displaystyle\Ket{3\pm}=\Ket{\Gamma_{6\pm}}=\ket{\pm\tfrac{1}{2}},$ (18)

where $w=0.95$[7] and the particle-number $n_{i\pm}$ of each KD
$|i\pm\rangle~{}(i=1,2,3)$ is also given by,

$\displaystyle
n_{1+}=n_{1-}=\frac{1}{Z}=\frac{1}{2(1+e^{-\Delta_{1}/T}+e^{-\Delta_{2}/T})},$
(19) $\displaystyle
n_{2+}=n_{2-}=\frac{e^{-\Delta_{1}/T}}{Z}=\frac{e^{-\Delta_{1}/T}}{2(1+e^{-\Delta_{1}/T}+e^{-\Delta_{2}/T})},$
(20) $\displaystyle
n_{3+}=n_{3-}=\frac{e^{-\Delta_{2}/T}}{Z}=\frac{e^{-\Delta_{2}/T}}{2(1+e^{-\Delta_{1}/T}+e^{-\Delta_{2}/T})}.$
(21)

The single-site susceptibility $\chi_{m_{1}m_{2}m_{3}m_{4}}^{(0)}$ becomes
finite only for certain sets of $m_{1},m_{2},m_{3},m_{4}$ as summarized in
Table 3, which results in the following form,

$\displaystyle\chi_{m_{1}m_{2}m_{3}m_{4}}^{(0)}=\delta_{m_{3}m_{1}}\delta_{m_{2}m_{4}}\Bigl{\\{}\delta_{11}^{m_{1}m_{2}}f_{1}+\delta_{22}^{m_{1}m_{2}}f_{2}+\delta_{33}^{m_{1}m_{2}}f_{3}$
$\displaystyle\qquad\qquad\qquad\qquad+\left(\delta_{12}^{m_{1}m_{2}}+\delta_{21}^{m_{1}m_{2}}\right)f_{4}+\left(\delta_{13}^{m_{1}m_{2}}+\delta_{31}^{m_{1}m_{2}}\right)f_{5}+\left(\delta_{23}^{m_{1}m_{2}}+\delta_{32}^{m_{1}m_{2}}\right)f_{6}\Bigr{\\}},$
(22)
$\displaystyle\quad\delta_{ij}^{m_{1}m_{2}}=\delta_{m_{1},i+}\delta_{m_{2},j+}+\delta_{m_{1},i-}\delta_{m_{2},j-}+\delta_{m_{1},i+}\delta_{m_{2},j-}+\delta_{m_{1},i-}\delta_{m_{2},j+}\quad(i,j=1,2,3),$
(23)

where $f_{1}$-$f_{6}$ are functions with different $T$-dependence given in
Table 3.

Figure 6 shows (a) the number of occupied electrons $n_{i}=n_{i+}+n_{i-}$ in
the three KDs and (b) the functions $f_{1}$-$f_{6}$ in the single-site
susceptibility as a function of $T$ under the present CEF. The occupied number
of the lowest KD $n_{1}$ (excited KDs $n_{2},n_{3}$) increases (decreases)
with decreasing $T$, and finally $n_{1}\rightarrow 1~{}(n_{2},n_{3}\rightarrow
0)$ for $T\lesssim 2~{}{\rm meV}$ as shown in Fig. 6(a). Corresponding to such
the $T$-dependence of $n_{i}$, the Curie term of the lowest KD $f_{1}$
increases in proportion to $1/T$ with decreasing $T$ as shown in Fig. 6(b).
The Van-Vleck terms between the lowest KD and the excited KDs, $f_{4},f_{5}$,
also increase with decreasing $T$, but they become constant at low
temperatures of $T\lesssim 2~{}{\rm meV}$. The $T,O_{\Gamma}$-dependent
differences in single-site multipole susceptibility shown in Fig. 2 of the
main text are due to the presence or absence of $f_{1}$ and $f_{4},f_{5}$ and
the size of their coefficients, for example, the $f_{1}$ contribution of
$H_{0}$ is larger than taht of $I$, resulting in the larger susceptibility of
$H_{0}$ than that of $I$.

Table 3:  Finite sets of the matrix elements of the single-site susceptibility
$\chi_{m_{1}m_{2}m_{3}m_{4}}^{(0)}$.

$\chi_{m_{1}m_{2}m_{3}m_{4}}^{(0)}$ | $m_{1},m_{2},m_{3},m_{4}$
---|---
$\displaystyle f_{1}\equiv\frac{1}{ZT}=\frac{n_{1\pm}}{T}$ |  | $1\pm,1\pm,1\pm,1\pm$
---
$1\pm,1\mp,1\pm,1\mp$
$\displaystyle f_{2}\equiv\frac{e^{-\Delta_{1}/T}}{ZT}=\frac{n_{2\pm}}{T}$ |  | $2\pm,2\pm,2\pm,2\pm$
---
$2\pm,2\mp,2\pm,2\mp$
$\displaystyle f_{3}\equiv\frac{e^{-\Delta_{2}/T}}{ZT}=\frac{n_{3\pm}}{T}$ |  | $3\pm,3\pm,3\pm,3\pm$
---
$3\pm,3\mp,3\pm,3\mp$
$\displaystyle f_{4}\equiv\frac{1-e^{-\Delta_{1}/T}}{Z\Delta_{1}}=\frac{n_{1\pm}-n_{2\pm}}{\Delta_{1}}$ |  | $1\pm,2\pm,1\pm,2\pm$
---
$1\pm,2\mp,1\pm,2\mp$
$2\pm,1\pm,2\pm,1\pm$
$2\pm,1\mp,2\pm,1\mp$
$\displaystyle f_{5}\equiv\frac{1-e^{-\Delta_{2}/T}}{Z\Delta_{2}}=\frac{n_{1\pm}-n_{3\pm}}{\Delta_{2}}$ |  | $1\pm,3\pm,1\pm,3\pm$
---
$1\pm,3\mp,1\pm,3\mp$
$3\pm,1\pm,3\pm,1\pm$
$3\pm,1\mp,3\pm,1\mp$
$\displaystyle f_{6}\equiv\frac{e^{-\Delta_{1}/T}-e^{-\Delta_{2}/T}}{Z(\Delta_{2}-\Delta_{1})}=\frac{n_{2\pm}-n_{3\pm}}{\Delta_{2}-\Delta_{1}}$ |  | $2\pm,3\pm,2\pm,3\pm$
---
$2\pm,3\mp,2\pm,3\mp$
$3\pm,2\pm,3\pm,2\pm$
$3\pm,2\mp,3\pm,2\mp$

##  4 Derivation of RKKY interactions

Here we derive the RKKY interaction
$K_{m_{1}m_{2}m_{3}m_{4}}^{\alpha,\beta}(\bm{q})$ of Eq. (3) in the main text.
First, we introduce the generalized periodic Anderson (GPAM) model[27] which
can be written by,

$\displaystyle\mathscr{H}_{\rm GPAM}$
$\displaystyle=\mathscr{H}_{c\textrm{-}c}+\mathscr{H}_{c\textrm{-}f}+\mathscr{H}_{f\textrm{-}f}$
(24)
$\displaystyle=\sum_{\bm{k}\ell\ell^{\prime}}h_{\ell\ell^{\prime}}^{cc}(\bm{k})c_{\bm{k}\ell}^{\dagger}c_{\bm{k}\ell^{\prime}}+\frac{1}{\sqrt{N}}\sum_{\bm{k}\ell}\sum_{\bm{i}\alpha
m}\left(e^{-i\bm{k}\cdot\bm{R_{i}}}V_{\bm{k}\ell\alpha
m}c_{\bm{k}\ell}^{\dagger}f_{\bm{i}\alpha
m}+h.c.\right)+\sum_{\bm{i}\alpha}\mathscr{H}^{f}_{\bm{i}\alpha},$ (25)

where the first term is the same in Eq. (10) and the second term corresponds
to the Fourier transform of the $f$-electron of the third term in Eq. (10),
where the sub-lattice degrees of freedom of Ce-$f$ electrons are explicitly
written by $\alpha={\rm Ce1,Ce2}$. The third term is the local $f$-electron
term including the onsite energy levels and Coulomb interactions for $f$
electrons at each Ce atoms, where we do not write them explicitly here.

Next, we derive the effective Hamiltonian restricted to $f^{1}$ states by a
second-order perturbation (SOP) w. r. t. $\mathscr{H}_{c\textrm{-}f}$ in
$\mathscr{H}_{\rm GPAM}$. The resulting Hamiltonian is the multi-orbital Kondo
lattice model which is explicitly given by,

$\displaystyle\mathscr{H}_{\rm MKL}$ $\displaystyle=\sum_{\bm{i}\alpha
m}\varepsilon_{m}^{f}f_{\bm{i}\alpha m}^{\dagger}f_{\bm{i}\alpha
m}+\sum_{\bm{k}\ell\ell^{\prime}}h_{\ell\ell^{\prime}}^{cc}(\bm{k})c_{\bm{k}\ell}^{\dagger}c_{\bm{k}\ell^{\prime}}+\sum_{\bm{i}\alpha
mm^{\prime}}\sum_{\bm{k}\bm{k^{\prime}}\ell\ell^{\prime}}J_{\bm{i}\alpha
mm^{\prime}}^{\bm{k}\ell,\bm{k^{\prime}}\ell^{\prime}}f_{\bm{i}\alpha
m}^{\dagger}f_{\bm{i}\alpha
m^{\prime}}c_{\bm{k}\ell}^{\dagger}c_{\bm{k^{\prime}}\ell^{\prime}},$ (26)

where $m$ represents 6 states of Eqs. (16)-(18) on each atom $\alpha$ in uc at
$\bm{R_{i}}$ and $\varepsilon_{m}^{f}$ is the CEF energy level of the $f^{1}$
state given by $O(10^{0-1}~{}{\rm meV})$ as written in the main text, while
the energy level of the $f^{0}$ state $\Delta_{0}$ is $O(10^{0}~{}{\rm eV})$
much larger than $\varepsilon_{m}^{f}$. Here we employ only the
$f^{0}$-intermediate state for the Kondo coupling $J_{\bm{i}\alpha
mm^{\prime}}^{\bm{k}\ell,\bm{k^{\prime}}\ell^{\prime}}$ in Eq. (26) and obtain
as follows,

$\displaystyle J_{\bm{i}\alpha
mm^{\prime}}^{\bm{k}\ell,\bm{k^{\prime}}\ell^{\prime}}=\frac{1}{N}V_{\bm{k}\ell\alpha
m}V_{\bm{k^{\prime}}\ell^{\prime}\alpha
m^{\prime}}^{*}\left(\frac{1}{\varepsilon_{\bm{k}\ell}^{c}+\Delta_{0}-\varepsilon_{m}^{f}}+\frac{1}{\varepsilon_{\bm{k^{\prime}}\ell^{\prime}}^{c}+\Delta_{0}-\varepsilon_{m^{\prime}}^{f}}\right)e^{-i(\bm{k}-\bm{k^{\prime}})\cdot\bm{R_{i}}},$
(27)

where the contribution of the $f^{2}$ intermediate states is treated as the
same as that of the $f^{0}$-intermediate state. In a study of the dynamical
mean field theory for CeB6[21], the $f$-electron susceptibility has been
calculated including the contribution of $f^{n}$-multiplet intermediate states
on 1-ion, where the contribution of $f^{2}$ states is found to be less
important for the main $\bm{q}$ dependence of the Stoner factor that
determines the critical point of the multipole order. Therefore the present
treatment taken only the $f^{0}$-contribution can give the proper
$\bm{q}$-dependence of the RKKY interactions and obtain the dominant
instability towards the multipole order. Moreover we fix the scattered
$c$-electron energies to the Fermi energy $(=0)$,
$\varepsilon_{\bm{k}\ell}^{c},~{}\varepsilon_{\bm{k^{\prime}}\ell^{\prime}}^{c}\rightarrow
0$. Then the Kondo coupling $J_{\bm{i}\alpha
mm^{\prime}}^{\bm{k}\ell,\bm{k^{\prime}}\ell^{\prime}}$ can be written by the
following form,

$\displaystyle J_{\bm{i}\alpha
mm^{\prime}}^{\bm{k}\ell,\bm{k^{\prime}}\ell^{\prime}}=\frac{1}{N}V_{\bm{k}\ell\alpha
m}V_{\bm{k^{\prime}}\ell^{\prime}\alpha
m^{\prime}}^{*}\left(\frac{1}{\Delta_{0}-\varepsilon_{m}^{f}}+\frac{1}{\Delta_{0}-\varepsilon_{m^{\prime}}^{f}}\right)e^{-i(\bm{k}-\bm{k^{\prime}})\cdot\bm{R_{i}}}.$
(28)

The RKKY Hamiltonian can be obtained from the SOP w. r. t. the Kondo coupling
term together with the thermal average for the $c$ states. The final form is
given by,

$\displaystyle\mathscr{H}_{\rm
RKKY}=-\sum_{\bm{i},\bm{\delta}}\sum_{m_{1}m_{2}}\sum_{m_{3}m_{4}}K_{m_{1}m_{2}m_{3}m_{4}}^{\alpha,\beta}(\bm{\delta})f_{\bm{i}\alpha
m_{1}}^{\dagger}f_{\bm{i}\alpha m_{2}}f_{\bm{i}+\bm{\delta}\beta
m_{4}}^{\dagger}f_{\bm{i}+\bm{\delta}\beta m_{3}},$ (29) $\displaystyle
K_{m_{1}m_{2}m_{3}m_{4}}^{\alpha,\beta}(\bm{\delta})=\frac{1}{N}\sum_{\bm{q}}K_{m_{1}m_{2}m_{3}m_{4}}^{\alpha,\beta}(\bm{q})~{}e^{i\bm{q}\cdot\bm{\delta}},$
(30)

where $K_{m_{1}m_{2}m_{3}m_{4}}^{\alpha,\beta}(\bm{\delta})$ is the RKKY
coupling between $f$ states $(\alpha m_{1}m_{2})$ at $\bm{R_{i}}$ and $(\beta
m_{3}m_{4})$ at $\bm{R_{i}}+\bm{\delta}$ with the primitive translation vector
$\bm{\delta}$. The key quantity
$K_{m_{1}m_{2}m_{3}m_{4}}^{\alpha,\beta}(\bm{q})$ is given by,

$\displaystyle
K_{m_{1}m_{2}m_{3}m_{4}}^{\alpha,\beta}(\bm{q})=C_{m_{1}m_{2}m_{3}m_{4}}\frac{1}{N}\sum_{nn^{\prime}\bm{k}}\mathcal{M}^{n\bm{k}}_{\beta
m_{3},\alpha m_{1}}\mathcal{M}^{n\bm{k}+\bm{q}}_{\alpha m_{2},\beta
m_{4}}\frac{f_{n^{\prime}\bm{k}+\bm{q}}-f_{n\bm{k}}}{\varepsilon_{n\bm{k}}-\varepsilon_{n^{\prime}\bm{k}+\bm{q}}},$
(31) $\displaystyle
C_{m_{1}m_{2}m_{3}m_{4}}=\left(\frac{1}{\Delta_{0}-\varepsilon_{m_{1}}^{f}}+\frac{1}{\Delta_{0}-\varepsilon_{m_{2}}^{f}}\right)\left(\frac{1}{\Delta_{0}-\varepsilon_{m_{3}}^{f}}+\frac{1}{\Delta_{0}-\varepsilon_{m_{4}}^{f}}\right),$
(32) $\displaystyle\mathcal{M}^{n\bm{k}}_{\alpha m,\beta
m^{\prime}}=\sum_{\ell\ell^{\prime}}V_{\bm{k}\ell\alpha
m}^{*}V_{\bm{k}\ell^{\prime}\beta
m^{\prime}}U_{n\bm{k}\ell}^{c*}U_{n\bm{k}\ell^{\prime}}^{c},$ (33)

where $C_{m_{1}m_{2}m_{3}m_{4}}$ is a energy denominator by the
$f^{1}$-$f^{0}$ process and $\mathcal{M}^{n\bm{k}}_{\alpha m,\beta
m^{\prime}}$ denotes the scattering matrix element between the $f$-electron
states $(\alpha m)$ and $(\beta m^{\prime})$ via $c$-$f$ mixing with the
$c$-electron band states with $(n\bm{k})$. Here we set the excitation energy
to $\Delta_{0}=2$ eV from the center of mass of the $f$ bands of LaCoSi shown
in Fig. 1(a) of the main text, which yields the energy denominator
$C_{m_{1}m_{2}m_{3}m_{4}}\simeq 4/\Delta_{0}^{2}=1~{}[{\rm eV^{-2}}]$ due to
$\Delta_{0}\gg\varepsilon_{m}^{f}$, resulting in the same form of
$K_{m_{1}m_{2}m_{3}m_{4}}^{\alpha,\beta}(\bm{q})$ as Eq. (3) in the main text.

##  5 Enhancement rates of the RPA susceptibilities

Figure 7 shows the enhancement rates of the $\bm{q}=\bm{0}$ RPA
susceptibilities for (a)-(d) uniform and (e)-(h) staggered components,
excluding the TR-even $\Gamma_{1}$ mulitpoles. It can be seen that the uniform
and staggered components exhibit similar behavior. The enhancement rates for
the quadrupoles and hexadecapoles have a small $T$-dependence and are hardly
enhanced as shown in Figs. 7(a),(e) and Figs. 7(b),(f). These behavior can be
explained by the $T$-dependence of the single-site susceptibilities in Fig.
2(a) and 2(b) of the main text, i.e., TR-even susceptibilities other than
$\Gamma_{1}$ do not include the lowest KD Curie term contribution and are at
most increased by the Van-Vleck term contribution. On the other hand, the
rates for TR-odd multipoles of
$(J_{z},~{}T_{z}^{\alpha},~{}D_{z}^{1\alpha,2\alpha})~{}[\Gamma_{2}]$ and
$(J_{x/y},~{}T_{x/y}^{\alpha},~{}D_{x/y}^{1\alpha,2\alpha})~{}[\Gamma_{5}]$
enhance two to three times as shown in Figs. 7(c),(d) and (g),(h), which is
attributed to the fact that the TR-odd susceptibilities shown in Figs. 2(c)
and (d) of the main text include a Curie term contribution.

Figure 7: (Color online) The enhancement rates of the $\bm{q}=\bm{0}$ RPA
susceptibility for (a)-(d) uniform and (e)-(h) staggered multipoles except for
the TR-even $\Gamma_{1}$ multipoles as a function of $T$ for (a),(e) monopole
and quadrupoles, (b),(f) hexadecapoles, (c),(g) dipoles and octupoles and
(d),(h) triakontadipoles.