# Multipolar ordering from dynamical mean field theory with application to
CeB6

Junya Otsuki  Research Institute for Interdisciplinary Science, Okayama
University, Okayama 700-8530, Japan  Kazuyoshi Yoshimi  Institute for Solid
State Physics, University of Tokyo, Chiba 277-8581, Japan  Hiroshi Shinaoka
Department of Physics, Saitama University, Saitama 338-8570, Japan  Harald O.
Jeschke  Research Institute for Interdisciplinary Science, Okayama University,
Okayama 700-8530, Japan

###### Abstract

Magnetic and multipolar ordering in $f$ electron systems takes place at low
temperatures of order 1–10 Kelvin. Combinations of first-principles with many-
body calculations for such low-energy properties of correlated materials are
challenging problems. We address multipolar ordering in $f$ electron systems
based on the dynamical mean-field theory (DMFT) combined with density
functional theory. We derive the momentum-dependent multipolar
susceptibilities and interactions by solving the Bethe-Salpeter (BS) equation
of the two-particle Green’s function. We apply the formalism to the
prototypical example of multipolar ordering CeB6, and demonstrate that the
experimental quadrupole transition is correctly reproduced. This first-
principles formalism based on DMFT and BS equation has applications which are
beyond the reach of the traditional RKKY formula. In particular, more
itinerant electron systems including $5f$, $4d$ and $5d$ electrons can be
addressed.

##  I Introduction

Recent progress in electronic structure calculations has been expanding the
target of first-principles calculations to strongly correlated materials, in
particular with the help of the DMFT [1, 2, 3]. The DMFT combined with density
functional theory (DFT+DMFT) has been employed extensively to investigate the
electronic structure of transition metal oxides [4, 5, 6, 7], of iron based
superconductors [8, 9, 10], of rare earth compounds [11, 12], and of many
other materials where electronic correlations in the form of local Coulomb
repulsion or Hund’s rule coupling are crucial. DFT+DMFT is often employed for
the theoretical support of photoemission experiments [13]. There is an ongoing
effort to account for the momentum dependence of the self-energy [14], and the
efficacy of new approaches is often demonstrated for the prime example of a
moderately correlated oxide, SrVO3. For example, DMFT has been combined with
the $GW$ method [15, 16], with vertex corrections [17, 18] or improved by the
dynamical cluster approximation [19].

The descriptions of phase transitions in correlated materials is one of the
objectives in DFT+DMFT investigations [20]. Typical examples include metallic
ferromagnetism in elemental Fe [21, 22, 23, 24, 25, 26, 27], isostructural
transitions in elemental Ce [28, 29, 30, 31, 32], and the orbital ordering in
Cu fluorides [33, 34]. Ordered magnetism has been studied within DFT+DMFT [35,
36, 37, 38, 39], but such investigations are less routine than the application
to paramagnetic phases. In particular, the determination of two-particle
susceptibilities for comparison with inelastic neutron scattering continues to
be difficult [40, 41, 42, 43].

Figure 1: A schematic $P$–$T$ phase diagram of $f$ electron materials and two
contrasting approaches to the multipolar ordering.

In this direction, multipolar ordering in rare-earth compounds is a
challenging subject because of the $f$ orbital degrees of freedom, the strong
spin-orbit coupling, and the small energy scale of order 10 K or even less. In
this paper, we address multipolar ordering based on the DFT+DMFT method.
Conventionally, the multipolar ordering is considered by RKKY interactions
where $f$ electrons are treated as localized. On the other hand, a nonmagnetic
Fermi liquid (heavy fermion) ground state can be addressed from the opposite
limit with itinerant $f$ electrons. Actually, the localized and itinerant $f$
states are continuously connected in a $T$-$P$ phase diagram (Fig. 1). DMFT
provides descriptions of the crossover between localized state and itinerant
state including the heavy-fermion state near the quantum critical point [44,
45, 46]. Therefore, formalism based on DMFT covers the whole region of the
phase diagram. This motivates us to establish descriptions of multipolar
ordering based on the DFT+DMFT framework for future wider applications such as
mixed-valent rare-earth compounds and actinide compounds with itinerant $5f$
electrons.

We focus on CeB6 as a prototypical material for multipolar ordering (for a
review, see Ref. [47]). The $4f^{1}$ configuration in trivalent Ce ions forms
$\Gamma_{8}$ quartet ground state of the total angular momentum $j=5/2$ under
the cubic crystalline electric field (CEF) [48]. Due to the four-fold
degeneracy of the CEF multiplet, CeB6 shows a rich phase diagram as shown in
Fig. 2. Two phase transitions have been observed at low temperatures [49].
Below $T_{\mathrm{Q}}=3.4$ K, CeB6 exhibits the antiferro-quadrupole (AFQ)
order of $\bm{q}=(1/2,1/2,1/2)$ (phase II) [50, 51, 52]. At
$T_{\mathrm{N}}=2.3$ K, the antiferro-magnetic (AFM) order takes place on top
of the AFQ order (phase III) [53]. Furthermore, another phase appears by
partially substituting Ce with La (phase IV) [54]. This phase is ascribed to
the antiferro-octupole order [55, 56, 57].

Figure 2: A schematic $T$–$H$ phase diagram of CeB6 ($\bm{H}\parallel[0,0,1]$)
[58, 54]. The lattice shows the configuration of the AFQ order in phase II,
where the color indicates the sign of the $xy$-type quadrupole.

Theoretically, the phase transitions in CeB6 are described by 15 multipoles in
$\Gamma_{8}$ quartet states [59]. Since $4f$ electrons are well localized in
CeB6 [60], a model consisting only of the local degrees of freedom is a good
starting point. Based on the Heisenberg-type model for multipoles, Shiina et
al. derived the $T$-$H$ phase diagram and observables related to the order
parameter [61, 62]. In this approach, the coupling constants are determined so
that the model yields results consistent with experiments.

Another important issue beyond phenomenology is the microscopic derivation of
the multipolar interactions to understand why the multipolar ordering takes
place in individual materials. Shiba et al. analyzed the RKKY formula in the
$\Gamma_{8}$ systems and deduced general properties [63]. Evaluation of the
multipolar interactions was performed later taking explicit electronic
structures into account for CeB6 [64, 65, 66, 67], NpO2 [68], and Pr compounds
[69].

Our approach based on the DFT+DMFT method includes the RKKY formula in the
strong-coupling limit [70]. In the following, we present an application of the
DFT+DMFT method to CeB6 and reproduce the experimental order parameter.
Through this calculation, we shall establish how to fix parameters such as the
Coulomb repulsion $U$ so that the method is applicable to other materials
without any tuning parameter.

##  II Single-particle properties

###  II.1 DFT calculations

The starting point of our study is a precise DFT calculation for CeB6. We
determine its fully relativistic electronic structure using the full potential
local orbital (FPLO) basis set [71]. The calculations are performed based on a
single crystal X-ray diffraction structure determined at $T=165$ K [72]. It
has a cubic space group $Pm\bar{3}m$ and features one Ce atom in the corner of
the cubic cell and a B6 octahedron in its center as shown in the inset of Fig.
3 (b). The DFT bands are shown in Fig. 3 (a) (full lines). We use symmetry
protected maximally projected Wannier functions [73, 74] to construct a 72
band model consisting of Ce $4f$ and $5d$ as well as B $2s$ and $2p$ orbitals;
the resulting tight binding representation of the bands is shown as dashed
lines in Fig. 3 (a). Figure 3 (b) shows the density of states (DOS). The Ce
$4f$ partial DOS is resolved by the total angular momentum $j$ contributions.
The energy splitting $\Delta_{\mathrm{SOC}}$ between $j=5/2$ and $j=7/2$ due
to the spin-orbit coupling (SOC) is $\Delta_{\mathrm{SOC}}=0.32$–0.34 eV,
which is comparable to the value in Ce ion. The CEF splitting between the
$\Gamma_{8}$ quartet ground state and the excited $\Gamma_{7}$ doublet is 11.3
meV ($\approx 131$ K), which is smaller by a factor of 4 compared to the
experimental value 540 K [48].

Figure 3: (a) Fully relativistic band structure of CeB6 (full lines), shown
together with the tight binding bands obtained from projective Wannier
funcions (dashed lines). High symmetry points of the cubic $Pm\bar{3}m$ space
group are $\Gamma=(0,0,0)$, ${\rm R}=(1/2,1/2,1/2)$, ${\rm M}=(1/2,1/2,0)$,
${\rm X}=(1/2,0,0)$. (b) Total and orbital-resolved density of states. The
$Pm\bar{3}m$ crystal structure of CeB6 is shown as inset.

###  II.2 DMFT calculations

Within the DFT+DMFT method, the single-particle Green’s function matrix
$\hat{G}(\bm{k},i\omega_{n})$ is given by [2]

$\displaystyle\hat{G}(\bm{k},i\omega_{n})=[(i\omega_{n}+\mu)\hat{I}-\hat{H}_{\mathrm{DFT}}(\bm{k})-\hat{\Sigma}_{\mathrm{loc}}(i\omega_{n})+\hat{\Sigma}_{\mathrm{DC}}]^{-1},$
(1)

where quantities with hat stand for $72\times 72$ matrices. $\omega_{n}$ is
the Matsubara frequency, $\mu$ is the chemical potential, and $\hat{I}$ is the
identity matrix. $\hat{H}_{\mathrm{DFT}}(\bm{k})$ denotes the non-interacting
tight-binding Hamiltonian constructed by the DFT calculations.
$\hat{\Sigma}_{\mathrm{loc}}(i\omega_{n})$ is the local self-energy in DMFT,
and $\hat{\Sigma}_{\mathrm{DC}}$ is the double-counting (DC) correction.

The local self-energy $\hat{\Sigma}_{\mathrm{loc}}(i\omega_{n})$ is computed
by solving the effective impurity problem for $f$ orbitals. We treat only
$j=5/2$ states, neglecting $j=7/2$ states ($j$-$j$ coupling scheme). The
influence of the neglected $j=7/2$ multiplet was analyzed in detail in Refs.
75, 76. We employ fully rotationally invariant Slater interactions, which are
specified by four parameters called Slater integrals $F_{0}$, $F_{2}$,
$F_{4}$, and $F_{6}$. We use the standard parameterization that relates the
Slater integrals to two intuitive parameters, the Coulomb repulsion $U$ and
the Hund’s coupling $J_{\mathrm{H}}$ [77].

We represent the DC self-energy in a simple form

$\displaystyle\hat{\Sigma}_{\mathrm{DC}}=-\Delta\epsilon_{f}\hat{N}_{f},$ (2)

where $\hat{N}_{f}$ is a projection onto $f$ orbitals. $\Delta\epsilon_{f}$
stands for the energy shift of the single-particle energy of the $f$
electrons. This expression is motivated by the fact that the local electronic
structure such as the SOC and the CEF are already taken into account by the
DFT and do not need to be changed. The $\hat{\Sigma}_{\mathrm{DC}}$ cancels
only the averaged Hartree shift in $\hat{\Sigma}_{\mathrm{loc}}(i\omega_{n})$
so that $4f$ electrons remain on the Ce atom.

We thus have 3 parameters to determine: $U$, $J_{\mathrm{H}}$, and
$\Delta\epsilon_{f}$. We determine these parameters so that the calculated
single-particle excitation spectrum agrees with the experimental and
theoretical PES and BIS spectra. Two energies characterize the spectrum:
$\Delta_{-}$ is the minimum excitation energy from the occupied $4f^{n}$
multiplet to the Fermi energy $E_{\mathrm{F}}$, and $\Delta_{+}$ is the
minimum excitation energy from $E_{\mathrm{F}}$ to empty $4f^{n+1}$
multiplets. Table 1 summarizes $\Delta_{-}$ and $\Delta_{+}$ in the
literature. From the PES experiment in CeB6, $\Delta_{-}$ has been reported to
be $\Delta_{-}=1.9$ eV. On the other hand, the BIS spectrum is not available
for CeB6. We instead refer the the data for elemental Ce. The experiment
reported $\Delta_{+}=3.46$ eV, while theoretical calculation yields
$\Delta_{+}=3.1$ eV. We adopt the theoretical value, since the theoretical
result for $\Delta_{-}$ in Ce agrees with the PES value for CeB6. From these
observations, we adopt $\Delta_{-}=1.9$ eV and $\Delta_{+}=3.1$ eV.

Table 1: Summary of the energy separation between occupied $4f^{1}$ state and Fermi level $\Delta_{-}$ and between occupied $4f^{1}$ state and lowest lying $4f^{2}$ state $\Delta_{+}$ in units of eV. material | $\Delta_{-}$ | $\Delta_{+}$
---|---|---
CeB6 (experiment) | 1.9 [78, 79] | –
Ce (experiment) | 0.27–1.92 [80] | 3.46 [80]
Ce (theory) | 1.9 [81] | 3.1 [82]

The DFT+DMFT cauculations were performed using open source software DCore [83]
implemented on TRIQS [84] and DFTTools [85] libraries. We solved the effective
impurity model by exact diagonalization. The exact diagonalization was
performed with pomerol [86]. Here, we neglected the effective hybridization
$\Delta(i\omega_{n})$. This corresponds to the Hubbard-I approximation, which
is often used for descriptions of properties of rare-earth compounds [87, 88,
89, 90]. In CeB6, the Fermi surface is quite similar to that of LaB6 [60],
indicating that the $4f$ electrons are well localized. Therefore, a reasonable
description by the Hubbard-I approximation is expected.

Figure 4: Eigenenergies of the many-body states of local $4f$ electrons. The
inset shows a zoom up of the low-energy multiplets.

Without $\Delta(i\omega_{n})$, energy levels of the many-body $4f$ states are
well defined. Figure 4 shows the eigenvalues of the $4f^{n}$ multiplets for
$U=6.2$ eV, $J_{\mathrm{H}}=0.8$ eV, and $\Delta\epsilon_{f}=-1.6$ eV. Our way
to determine these values is discussed in detail in Appendix A. The
corresponding Slater integrals are $F_{0}=6.2$ eV, $F_{2}=9.54$ eV,
$F_{4}=6.37$ eV, and $F_{6}=4.71$ eV. There are $2^{6}=64$ levels in this
figure. The CEF splitting is not visible on this scale, and we discuss only
the multiplets labeled by the total angular momentum $J$. The energy of the
$j=5/2$ multiplet of the $4f^{1}$ configuration is $E_{1}=-1.9$ eV
($\Delta_{-}=-E_{1}$). The $4f^{2}$ configuration has three multiplets, $J=4$,
2, 0. The Hund’s ground state $J=4$ is located at $E_{2}=+1.2$ eV, thus
$\Delta_{+}=E_{2}-E_{1}=3.1$ eV.

Figure 5: (a) The single-particle excitation spectrum $A(\bm{k},\omega)$
computed within DMFT. Only $j=5/2$ states are retained. (b) The
$\bm{k}$-integrated spectrum $A(\omega)$ at a temperature $T=0.01$ eV.

The final results for the single-particle excitation spectrum
$A(\bm{k},\omega)$ and the $\bm{k}$-integrated spectrum $A(\omega)$ are shown
in Fig. 5. The $4f$ peaks near $E_{\mathrm{F}}$ in the DFT result [Fig. 3 (b)]
are moved away from $E_{\mathrm{F}}$ in the presence of interactions. The peak
at $-1.9$ eV corresponds to $j=5/2$ multiplet of the $4f^{1}$ configuration,
and the peak at 3.1 eV corresponds to the $J=4$ multiplet of the $4f^{2}$
configuration. The resultant energy dispersion near $E_{\mathrm{F}}$ is
similar to LaB6 [60, 91]. A parabolic band around the X point forms a Fermi
surface that are connected on the $\Gamma$–M line. The occupied part of our
spectral function, Fig. 5 (a), agrees very well with photoemission experiment
[92] and is therefore a good starting point for the determination of two-
particle quantities.

##  III Two-particle properties

###  III.1 Multipolar susceptibility

Figure 6: The $\bm{q}$-dependence of 36 eigenvalues of the susceptibility
matrix $\hat{\chi}(\bm{q})$ in $j=5/2$ states. $T=0.01$ eV.

The static susceptibility describing fluctuations within $j=5/2$ states is
given by

$\displaystyle\chi_{m_{1}m_{2},m_{3}m_{4}}(\bm{q})=\int_{0}^{\beta}d\tau\langle
O_{m_{1}m_{2}}(\bm{q},\tau)O_{m_{4}m_{3}}(-\bm{q})\rangle.$ (3)

Here, the argument $\tau$ stands for the Heisenberg operator. The operator
$O_{mm^{\prime}}(\bm{q})$ is the Fourier transform of the local density
operator defined by

$\displaystyle O_{mm^{\prime}}(i)=f_{im}^{{\dagger}}f_{im^{\prime}},$ (4)

where $f_{im}^{{\dagger}}$ and $f_{im}$ are the creation and annihilation
operator for $f$ electrons, respectively. The subscript $m$ stands for the
eigenvalues of $j_{z}$, namely, $m=-5/2$, $-3/2$, $\cdots$, $+5/2$. Therefore,
there are $6^{2}=36$ components in the density operators, and $36^{2}=1{,}296$
components in the susceptibility. Regarding $O_{mm^{\prime}}$ as a vector,
$\chi_{m_{1}m_{2},m_{3}m_{4}}$ can be regarded as a $(36\times 36)$ matrix,
which we denote by $\hat{\chi}(\bm{q})$.

In DMFT, $\hat{\chi}(\bm{q})$ can be computed by solving the BS equation of
the two-particle Green’s function [93, 1]. The two-particle Green’s function
depends on two fermionic Matsubara frequencies in addition to four $m$
indices. We introduced a cutoff $\omega_{\mathrm{max}}$ for the frequency and
solved the matrix equation of size $(36N_{\omega}\times 36N_{\omega})$, where
$N_{\omega}$ is the number of Matsubara frequencies below
$\omega_{\mathrm{max}}$. All figures presented hereafter are results computed
with $N_{\omega}=160$ at $T=0.01$ eV ($\omega_{\mathrm{max}}\approx 5.1$ eV),
unless otherwise specified. For quantitative discussions, we extrapolate
results to $\omega_{\mathrm{max}}\to\infty$ using results computed with up to
$N_{\omega}=400$.

Figure 6 shows the eigenvalues of the susceptibility matrix. There are 36
eigenvalues, which are classified into 5 groups. The leading group contains
$15=4^{2}-1$ eigenvalues, which correspond to multipolar fluctuations within
$\Gamma_{8}$ quartet states. The second group contains 16 eigenvalues, which
arise due to mixing between $\Gamma_{8}$ and $\Gamma_{7}$ states (two times
$4\times 2$). The third largest is the single fluctuation mode with A1g
symmetry (hexadecapole). The fourth largest (second smallest) eigenvalues
group contains $3=2^{2}-1$ fluctuation modes, which correspond to the spin
fluctuations within $\Gamma_{7}$ doublet states. The lowest (almost zero)
single-mode fluctuation is the charge fluctuation.

The eigenvectors of the susceptibility matrix yield multipoles that take the
crystal symmetry into account. Instead, we can take linear combinations of
$O_{mm^{\prime}}$ using pre-computed coefficients
$C^{(\gamma)}_{mm^{\prime}}$, namely,

$\displaystyle O_{\gamma}=C^{(\gamma)}_{mm^{\prime}}O_{mm^{\prime}},$ (5)

where the subscript $\gamma$ distinguishes the irreducible representations. In
point group Oh, 15 operators in the $\Gamma_{8}$ quartet system (excluding the
charge) are classified into 6 kinds of multipoles according to the rank and
irreducible representations [62, 94, 47]. Table 2 summarizes the basis
function of these multipoles. The corresponding expressions for
$C^{(\gamma)}_{mm^{\prime}}$ can be constructed from the basis set [62, 94,
47]. We adopted a normalization
$\sum_{mm^{\prime}}|C^{(\gamma)}_{mm^{\prime}}|^{2}=1$. The multipolar
susceptibilities $\chi_{\gamma}(\bm{q})$ can then be evaluated by taking
linear combinations of $\chi_{m_{1}m_{2},m_{3}m_{4}}(\bm{q})$ according to

$\displaystyle\chi_{\gamma}(\bm{q})=\sum_{m_{1}m_{2}m_{3}m_{4}}C_{m_{1}m_{2}}^{(\gamma)}\chi_{m_{1}m_{2},m_{3}m_{4}}(\bm{q})C_{m_{3}m_{4}}^{(\gamma)\ast}.$
(6)

Table 2: The basis functions of 15 multipoles in the $\Gamma_{8}$ quartet system classified by the rank and irreducible representations (irreps). Extracted from Ref. [47]. | Irrep | Basis
---|---|---
Dipole | $\Gamma_{4u}$ | $\\{x,\ y,\ z\\}$
Quadrupole | $\Gamma_{3g}$ | $\\{3z^{2}-r^{2},\ \sqrt{3}(x^{2}-y^{2})\\}$
| $\Gamma_{5g}$ | $\\{xy,\ yz,\ zx\\}$
Octupole | $\Gamma_{2u}$ | $xyz$
| $\Gamma_{4u}$ | $\\{x(5x^{2}-3r^{2}),\ y(5y^{2}-3r^{2}),\ z(5z^{2}-3r^{2})\\}$
| $\Gamma_{5u}$ | $\\{x(y^{2}-z^{2}),\ y(z^{2}-x^{2}),\ z(x^{2}-y^{2})\\}$
Figure 7: The $\bm{q}$ dependence of 15 multipolar susceptibilities
$\chi_{\gamma}(\bm{q})$ within $\Gamma_{8}$ states.

Figure 7 shows 15 susceptibilities $\chi_{\gamma}(\bm{q})$ within the
$\Gamma_{8}$ states. The leading fluctuations are the quadrupoles of
$\Gamma_{5g}$ symmetry at the R point
[$\bm{q}=(1/2,1/2,1/2)\equiv\bm{q}_{\mathrm{R}}$]. This is consistent with the
AFQ order of phase II in CeB6. The second largest fluctuation is the octupole
fluctuation of $\Gamma_{2u}$ type, which is almost degenerate with the leading
quadrupolar fluctuation. We will discuss this degeneracy later in the context
of the multipolar interactions.

In order to determine the transition temperature, we need to decrease $T$ to
detect a divergence of $\chi_{\gamma}(\bm{q})$ (see for example Ref. [23]).
However, it is practically difficult to perform the computations down to the
expected transition temperature of order 1–10 K because the required number of
Matsubara frequency, $N_{\omega}$, grows in inversely proportional to $T$.
Hence, we adopt another strategy, namely, of extracting the effective
intersite interactions.

###  III.2 Multipolar interactions

If the $4f$ electrons are well localized, we expect that a description using
an effective Heisenberg-type model with multipolar degrees of freedom is a
reasonable approximation. In this case, $\hat{\chi}(\bm{q})$ computed by
solving the BS equation follows the form
$\hat{\chi}(\bm{q})\approx[\hat{\chi}_{\mathrm{loc}}^{-1}-\hat{I}(\bm{q})]^{-1}$,
since DMFT treats intersite correlations at a mean-field level. Assuming this
expression, we _define_ the effective intersite interaction $\hat{I}(\bm{q})$
by

$\displaystyle\hat{I}(\bm{q})\equiv\hat{\chi}_{\mathrm{loc}}^{-1}-\hat{\chi}(\bm{q})^{-1}.$
(7)

Once $\hat{I}(\bm{q})$ is obtained, we can extrapolate $\hat{\chi}(\bm{q})$ to
lower temperatures because the dominant temperature dependence arises from
$\hat{\chi}_{\mathrm{loc}}\propto 1/T$ in localized $4f$ electron systems. We
note that in the actual evaluation of $\hat{I}(\bm{q})$, the charge
fluctuation has to be eliminated to avoid numerical instability in the matrix
inversion (see Appendix B).

Figure 8: The $\bm{q}$ dependence of the multipolar interactions
$I_{\gamma}(\bm{q})$ within $\Gamma_{8}$ states. Six kinds of multipoles are
plotted separately. From (a) to (f) correspond to top to bottom in Table 2.

Figure 8 shows $I_{\gamma}(\bm{q})$ plotted separately for each irreducible
representation. Here, positive (negative) values enhance (suppress) the
fluctuations. The largest interaction is the quadrupole of $\Gamma_{5g}$ type
at the R point, which induces the AFQ order in phase II. We express this
interaction by $I_{\mathrm{Q}}(\bm{q}_{\mathrm{R}})$. Figure 8 indicates that
the interaction strength of $I_{\mathrm{Q}}(\bm{q}_{\mathrm{R}})$ is about 1
meV. An extrapolation to $\omega_{\mathrm{max}}\to\infty$ yields
$I_{\mathrm{Q}}(\bm{q}_{\mathrm{R}})\approx 1.9$ meV (see Appendix C for
details). The transition temperature can be determined by the condition

$\displaystyle\chi_{\mathrm{loc}}I_{\mathrm{Q}}(\bm{q}_{\mathrm{R}})=1.$ (8)

Within the Hubbard-I approximation, the local susceptibility
$\chi_{\mathrm{loc}}$ follows the Curie law $\chi_{\mathrm{loc}}=1/4T$, where
we used $\sum_{mm^{\prime}}|C^{(\gamma)}_{mm^{\prime}}|^{2}=1$ and the factor
$4$ is due to the dimension of the $\Gamma_{8}$ states. Then, the transition
temperature $T_{\mathrm{Q}}$ is estimated to be
$T_{\mathrm{Q}}=I_{\mathrm{Q}}(\bm{q}_{\mathrm{R}})/4\approx
0.48\,\textrm{meV}\approx 5.6$ K. This is consistent with the experimental
value $T_{\mathrm{Q}}=3.4$ K.

The second largest interaction is the octupole of $\Gamma_{2u}$ type. This
result is consistent with the general conclusion in Ref. 63: “If the AFQ
interaction of $\Gamma_{5g}$ type is strong, the antiferro-octupolar
interaction of $\Gamma_{2u}$ type is equally strong”. The $\Gamma_{2u}$-type
antiferro-octupolar interaction plays a relevant role in stabilizing the AFQ
ordering in phase II under magnetic field [62, 47].

The order parameter in phase IV is reported as the octupole of $\Gamma_{5u}$
type with $\bm{q}=\bm{q}_{\mathrm{R}}$ [47]. Our result shows that the
octupolar interaction of $\Gamma_{5u}$ type is largest at the X point and
second largest at the R point. Hence, the interaction for phase IV is not
relevant in pure CeB6 without La substitution.

##  IV Doping dependence

In order to clarify which energy band plays a major role in the experimental
quadrupolar ordering, we artificially shift the chemical potential and
investigate how the effective interactions $I_{\gamma}(\bm{q})$ vary. Such a
“numerical experiment” is one of the advantages of theoretical calculations.
We shifted the chemical potential by up to $\pm 1$ eV in a rigid band
treatment. The shift of the chemical potential changes the electron number per
unit cell by $-0.6$ for $\Delta\mu=-1$ eV and $+0.9$ for $\Delta\mu=+1$ eV,
while the $4f$ electron number is unchanged. Therefore, electron (hole) doping
per boron is 0.15 (0.10) at $\Delta\mu=+1$ eV ($-1$ eV). The Fermi surface are
changed as follows: Under the hole doping, the open Fermi surface centered at
X point are closed (the dispersion on $\Gamma$–M line in Fig. 5 leaves away
from $E_{\mathrm{F}}$). On the other hand, the electron doping generates small
electron pockets on $\Gamma$–M line, and hence another contribution to the
susceptibility is expected.

Figure 9 shows the variation of the effective interactions
$I_{\gamma}(\bm{q})$ as a function of $\Delta\mu$. We plotted only relevant
interactions: the $\Gamma_{5g}$-quadrupole and the $\Gamma_{2u}$-octupole at R
point and $\Gamma$ point, and the $\Gamma_{3g}$-quadrupole and the
$\Gamma_{4u}$-octupole at the M point. The leading interaction
$I_{\mathrm{Q}}(\bm{q}_{\mathrm{R}})$ for $\Delta\mu=0$ is enhanced by the
hole doping ($\Delta\mu<0$). Therefore, we conclude that the Fermi surface
around the X point causes the experimental quadrupolar ordering. Upon further
hole doping, the octupolar interaction of $\Gamma_{4u}$ type becomes dominant.

On the other hand, $I_{\mathrm{Q}}(\bm{q}_{\mathrm{R}})$ is suppressed by
electron doping ($\Delta\mu>0$), which leads to a reduction of the transition
temperature of phase II. Instead, the same multipole with a different
$\bm{q}$-vector, namely, the $\Gamma_{5g}$-type ferro-quadrupolar interaction
$I_{\mathrm{Q}}(\bm{0})$ is enhanced. Upon further electron doping, the
quadrupolar interaction of $\Gamma_{3g}$ type at the M point takes the place
of the leading interaction. Our theory thus predicts that the the quadrupolar
order in phase II is robust against hole doping but is replaced by other order
parameters upon electron doping.

Figure 9: Variations of relevant effective interactions $I_{\gamma}(\bm{q})$
against the shift of the chemical potential, $\Delta\mu$. The symbols indicate
the $\bm{q}$-vector: closed circles for the R point, open circles for the
$\Gamma$ point, and closed triangles for the M point.

##  V Summary

We demonstrated a microscopic derivation of multipolar ordering using the
DFT+DMFT method. There are three parameters in these calculations: $U$,
$J_{\mathrm{H}}$, and $\Delta\epsilon_{f}$. We showed a recipe to fix them:
$\epsilon_{f}$ is determined from the PES spectrum, and $U-J_{\mathrm{H}}$ can
be determined from BIS spectrum. By fixing $J_{\mathrm{H}}$ to a typical
value, we can thus determine these parameters in an unbiased manner.

We evaluated the $\bm{q}$-dependent multipolar susceptibilities and
interactions in CeB6 by solving the BS equation. Although direct observation
of the second-order phase transition is practically difficult at present, we
can estimate the transition temperature by extrapolating the susceptibility to
low temperatures. Our results yield a good estimate.

Comparing our approach to the more traditional combination of multiorbital
Kondo lattice model and RKKY formula, there are a few advantages to the new
method: (i) We base our calculations on the full DFT+DMFT spectral function so
that the validity of our calculation can already be checked at the single
particle level by comparing to photoemission spectroscopy experiment. (ii) Our
approach is free from highly compound specific parameters that need to be
determined from experiment; with only two local interaction parameters $U$ and
$J_{\rm H}$ and a binding energy of the single $f$ electron, the method can be
applied to numerous Ce compounds. (iii) Our calculation procedure based on the
DFT+DMFT method is not limited to highly localized $f$ electron systems but
can be also applied to itinerant $f$ electron systems. The only difference to
the calculations in this paper is the impurity solver: We have to use more
sophisticated impurity solvers such as the continuous-time quantum Monte Carlo
method [95] to incorporate the formation of heavy-fermion states. Interesting
applications include a hidden order in URu2Si2 and scalar order in PrFe4P12
[47]. These are two of the many possible future uses of our method.

###### Acknowledgements.

This work was supported by JSPS KAKENHI grants No. 18H04301 (J-Physics), No.
20K20522, No. 21H01003, and No. 21H01041.

##  Appendix A Determination of parameters $U$, $J_{\mathrm{H}}$, and
$\Delta\epsilon_{f}$

We present the procedure of determining three parameters $U$,
$J_{\mathrm{H}}$, and $\Delta\epsilon_{f}$ so that $\Delta_{-}$ and
$\Delta_{+}$ agree the target values. We first vary $\Delta\epsilon_{f}$ with
$U$ and $J_{\mathrm{H}}$ fixed at 6.0 eV and 0.8 eV. Figure 10 shows
$\Delta_{-}$ and $\Delta_{+}$ as a function of $\Delta\epsilon_{f}$. The lower
graph for $\Delta_{-}$ is plotted with an inverted $y$-axis because of the
definition $\Delta_{-}=-E_{1}$. We find that the target value $\Delta_{-}=1.9$
eV is obtained at $\Delta\epsilon_{f}=-1.6$ eV. For information, when we
estimate $\Delta\epsilon_{f}$ from the Hartree energy, we obtained
$\Delta\epsilon_{f}\approx-5.3$ eV, which yields unphysically deep $4f$ level.

Figure 10: Calculated $\Delta_{-}$ and $\Delta_{+}$ as a function of
$\Delta\epsilon_{f}$ for $U=6.0$ eV and $J_{\mathrm{H}}=0.8$ eV. The
horizontal dashed lines indicate the target values, $\Delta_{-}=1.9$ eV and
$\Delta_{+}=3.1$ eV. The vertical dashed line indicates the value of
$\Delta\epsilon_{f}$ estimated by the Hartree energy.

Next, we turn our attention to $\Delta_{+}$. The sum $\Delta_{-}+\Delta_{+}$
(the difference between $4f^{1}$ and $4f^{2}$ peaks in $A(\omega)$) depends on
$U-J_{\mathrm{H}}$, because $\Delta_{+}$ is the excitation energy from
$E_{\mathrm{F}}$ to the Hund’ ground state of the $4f^{2}$ configuration.
Hence, the choice of $U$ and $J_{\mathrm{H}}$ themselves are arbitrary.
Following Ref. [89], we fix $J_{\mathrm{H}}=0.8$ eV and vary $U$. Then, we
found that $U=6.2$ eV gives the target value $\Delta_{+}=3.1$ eV.

##  Appendix B Calculation of the effective interaction

Direct calculation of Eq. (7) is unstable because the charge fluctuation is
tiny (or zero within the computer accuracy) and the inversion of the matrices
$\hat{\chi}_{\mathrm{loc}}$ and $\hat{\chi}_{\mathrm{loc}}$ diverges.
Actually, this divergence is canceled out to yield finite result. We present
here how to avoid this instability.

We decompose $\chi_{\mathrm{loc}}$ using the singular value decomposition
(SVD). Since $\chi_{\mathrm{loc}}$ is hermitian, SVD is reduced to
diagonalization

$\displaystyle\hat{\chi}_{\mathrm{loc}}=USU^{{\dagger}},$ (9)

where $U$ is a unitary matrix, and $S$ is a diagonal matrix whose diagonal
elements (singular values) are positive and aligned in descending order. It is
clear from the above expression that small singular values make the matrix
inversion unstable. Suppose that the smallest singular value is zero, then we
eliminate it and the corresponding basis from $U$. In this singular space, $S$
and $U$ are replaced, respectively, with an $[(M-1)\times(M-1)]$ diagonal
matrix $\tilde{S}$ and $[M\times(M-1)]$ matrix $\tilde{U}$. Eq. (7) is then
expressed as

$\displaystyle\hat{I}(\bm{q})\simeq\tilde{U}[\tilde{S}^{-1}-(\tilde{U}^{{\dagger}}\hat{\chi}(\bm{q})\tilde{U})^{-1}]\tilde{U}^{{\dagger}}.$
(10)

This expression does not suffer from numerical instability.

##  Appendix C Extrapolation of the frequency cutoff

Figure 11: Four leading effective interactions $I_{\gamma}(\bm{q})$ at
$\bm{q}=\bm{q}_{\mathrm{R}}$ as a function of $1/\omega_{\mathrm{max}}$. The
dashed lines show the result of linear fitting.

We solved the BS equation for $\chi({\bm{q}})$ by introducing the cutoff
$\omega_{\mathrm{max}}$ for the fermionic Matsubara frequency. In this
Appendix, we present extrapolation to $\omega_{\mathrm{max}}\to\infty$. Figure
11 shows the effective interaction $I_{\gamma}(\bm{q}_{\mathrm{R}})$ as a
function of $1/\omega_{\mathrm{max}}$. The data are well approximated by a
line as indicated by the dashed lines. Thus, extrapolation to
$1/\omega_{\mathrm{max}}\to 0$ yields $I_{\gamma}(\bm{q}_{\mathrm{R}})\approx
1.93$ meV for the $\Gamma_{5g}$-type quadrupolar interaction.

More elaborate methods for solving the BS equation have been proposed [40, 96,
97, 98, 99]. These approaches make the convergence faster or may allow us to
obtain converged results without extrapolation.

## References

  * Georges  _et al._ [1996] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68, 13 (1996).
  * Kotliar  _et al._ [2006] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Electronic structure calculations with dynamical mean-field theory, Rev. Mod. Phys. 78, 865 (2006).
  * Imada and Miyake [2010] M. Imada and T. Miyake, Electronic Structure Calculation by First Principles for Strongly Correlated Electron Systems, J. Phys. Soc. Jpn. 79, 112001 (2010).
  * Held [2007] K. Held, Electronic structure calculations using dynamical mean field theory, Adv. Phys. 56, 829 (2007).
  * Gorelov  _et al._ [2010] E. Gorelov, M. Karolak, T. O. Wehling, F. Lechermann, A. I. Lichtenstein, and E. Pavarini, Nature of the Mott transition in Ca2RuO4, Phys. Rev. Lett. 104, 226401 (2010).
  * Backes  _et al._ [2016] S. Backes, T. C. Rödel, F. Fortuna, E. Frantzeskakis, P. Le Fèvre, F. Bertran, M. Kobayashi, R. Yukawa, T. Mitsuhashi, M. Kitamura, K. Horiba, H. Kumigashira, R. Saint-Martin, A. Fouchet, B. Berini, Y. Dumont, A. J. Kim, F. Lechermann, H. O. Jeschke, M. J. Rozenberg, R. Valentí, and A. F. Santander-Syro, Hubbard band versus oxygen vacancy states in the correlated electron metal SrVO3, Phys. Rev. B 94, 241110 (2016).
  * Ahn  _et al._ [2021] C. Ahn, A. Cavalleri, A. Georges, S. Ismail-Beigi, A. J. Millis, and J.-M. Triscone, Designing and controlling the properties of transition metal oxide quantum materials, Nat. Mater. 20, 1462 (2021).
  * Aichhorn  _et al._ [2009] M. Aichhorn, L. Pourovskii, V. Vildosola, M. Ferrero, O. Parcollet, T. Miyake, A. Georges, and S. Biermann, Dynamical mean-field theory within an augmented plane-wave framework: Assessing electronic correlations in the iron pnictide LaFeAsO, Phys. Rev. B 80, 085101 (2009).
  * Yin  _et al._ [2011] Z. P. Yin, K. Haule, and G. Kotliar, Magnetism and charge dynamics in iron pnictides, Nat. Phys. 7, 294 (2011).
  * Ferber  _et al._ [2012] J. Ferber, H. O. Jeschke, and R. Valentí, Fermi surface topology of LaFePO and LiFeP, Phys. Rev. Lett. 109, 236403 (2012).
  * Haule  _et al._ [2010] K. Haule, C.-H. Yee, and K. Kim, Dynamical mean-field theory within the full-potential methods: Electronic structure of CeIrIn5, CeCoIn5, and CeRhIn5, Phys. Rev. B 81, 195107 (2010).
  * Lu and Huang [2016] H. Lu and L. Huang, Pressure-driven $4f$ localized-itinerant crossover in heavy-fermion compound ${\mathrm{CeIn}}_{3}$: A first-principles many-body perspective, Phys. Rev. B 94, 075132 (2016).
  * Ebert  _et al._ [2011] H. Ebert, D. Ködderitzsch, and J. Minár, Calculating condensed matter properties using the KKR-Green's function method—recent developments and applications, Rep. Prog. Phys. 74, 096501 (2011).
  * Rohringer  _et al._ [2018] G. Rohringer, H. Hafermann, A. Toschi, A. A. Katanin, A. E. Antipov, M. I. Katsnelson, A. I. Lichtenstein, A. N. Rubtsov, and K. Held, Diagrammatic routes to nonlocal correlations beyond dynamical mean field theory, Rev. Mod. Phys. 90, 025003 (2018).
  * Sakuma  _et al._ [2013] R. Sakuma, P. Werner, and F. Aryasetiawan, Electronic structure of SrVO3 within $GW$+DMFT, Phys. Rev. B 88, 235110 (2013).
  * Tomczak  _et al._ [2014] J. M. Tomczak, M. Casula, T. Miyake, and S. Biermann, Asymmetry in band widening and quasiparticle lifetimes in SrVO3: Competition between screened exchange and local correlations from combined $GW$ and dynamical mean-field theory $GW+\mathrm{DMFT}$, Phys. Rev. B 90, 165138 (2014).
  * Meng  _et al._ [2014] Z. Y. Meng, Y. B. Kim, and H.-Y. Kee, Odd-Parity Triplet Superconducting Phase in Multiorbital Materials with a Strong Spin-Orbit Coupling: Application to Doped ${\mathrm{Sr}}_{2}{\mathrm{IrO}}_{4}$, Phys. Rev. Lett. 113, 177003 (2014).
  * Galler  _et al._ [2017] A. Galler, P. Thunström, P. Gunacker, J. M. Tomczak, and K. Held, Ab initio dynamical vertex approximation, Phys. Rev. B 95, 115107 (2017).
  * Lee  _et al._ [2012] H. Lee, K. Foyevtsova, J. Ferber, M. Aichhorn, H. O. Jeschke, and R. Valentí, Dynamical cluster approximation within an augmented plane wave framework: Spectral properties of srvo3, Phys. Rev. B 85, 165103 (2012).
  * Kuneš  _et al._ [2017] J. Kuneš, I. Leonov, P. Augustinský, V. Křápek, M. Kollar, and D. Vollhardt, LDA+DMFT approach to ordering phenomena and the structural stability of correlated materials, Eur. Phys. J. Special Topics 226, 2641 (2017).
  * Lichtenstein  _et al._ [1987] A. Lichtenstein, M. Katsnelson, V. Antropov, and V. Gubanov, Local spin density functional approach to the theory of exchange interactions in ferromagnetic metals and alloys, J. Magn. Magn. Mater. 67, 65 (1987).
  * Halilov  _et al._ [1998] S. V. Halilov, H. Eschrig, A. Y. Perlov, and P. M. Oppeneer, Adiabatic spin dynamics from spin-density-functional theory: Application to Fe, Co, and Ni, Phys. Rev. B 58, 293 (1998).
  * Lichtenstein  _et al._ [2001] A. I. Lichtenstein, M. I. Katsnelson, and G. Kotliar, Finite-Temperature Magnetism of Transition Metals: An ab initio Dynamical Mean-Field Theory, Phys. Rev. Lett. 87, 067205 (2001).
  * Leonov  _et al._ [2011] I. Leonov, A. I. Poteryaev, V. I. Anisimov, and D. Vollhardt, Electronic Correlations at the $\alpha\mathrm{\text{-}}\gamma$ Structural Phase Transition in Paramagnetic Iron, Phys. Rev. Lett. 106, 106405 (2011).
  * Belozerov  _et al._ [2013] A. S. Belozerov, I. Leonov, and V. I. Anisimov, Magnetism of iron and nickel from rotationally invariant hirsch-fye quantum monte carlo calculations, Phys. Rev. B 87, 125138 (2013).
  * Szilva  _et al._ [2013] A. Szilva, M. Costa, A. Bergman, L. Szunyogh, L. Nordström, and O. Eriksson, Interatomic exchange interactions for finite-temperature magnetism and nonequilibrium spin dynamics, Phys. Rev. Lett. 111, 127204 (2013).
  * Nomoto  _et al._ [2020] T. Nomoto, T. Koretsune, and R. Arita, Local force method for the ab initio tight-binding model: Effect of spin-dependent hopping on exchange interactions, Phys. Rev. B 102, 014444 (2020).
  * Pickett  _et al._ [1981] W. E. Pickett, A. J. Freeman, and D. D. Koelling, Local-density-functional approach to the isostructual $\gamma$-$\alpha$ transition in cerium using the self-consistent linearized-augmented-plane-wave method, Phys. Rev. B 23, 1266 (1981).
  * Held  _et al._ [2001] K. Held, A. K. McMahan, and R. T. Scalettar, Cerium Volume Collapse: Results from the Merger of Dynamical Mean-Field Theory and Local Density Approximation, Phys. Rev. Lett. 87, 276404 (2001).
  * McMahan  _et al._ [2003] A. K. McMahan, K. Held, and R. T. Scalettar, Thermodynamic and spectral properties of compressed Ce calculated using a combined local-density approximation and dynamical mean-field theory, Phys. Rev. B 67, 075108 (2003).
  * Sakai  _et al._ [2005] O. Sakai, Y. Shimizu, and Y. Kaneta, Band Calculation for Ce-Compounds on the basis of Dynamical Mean Field Theory, J. Phys. Soc. Jpn. 74, 2517 (2005).
  * Lanatà  _et al._ [2013] N. Lanatà, Y.-X. Yao, C.-Z. Wang, K.-M. Ho, J. Schmalian, K. Haule, and G. Kotliar, $\gamma\mathrm{\text{$-$}}\alpha$ isostructural transition in cerium, Phys. Rev. Lett. 111, 196801 (2013).
  * Pavarini  _et al._ [2008] E. Pavarini, E. Koch, and A. I. Lichtenstein, Mechanism for Orbital Ordering in ${\mathrm{KCuF}}_{3}$, Phys. Rev. Lett. 101, 266405 (2008).
  * Musshoff  _et al._ [2019] J. Musshoff, G. Zhang, E. Koch, and E. Pavarini, Linear-response description of superexchange-driven orbital ordering in $\mathrm{K}{}_{2}\mathrm{CuF}{}_{4}$, Phys. Rev. B 100, 045116 (2019).
  * Matsumoto  _et al._ [2009] M. Matsumoto, M. J. Han, J. Otsuki, and S. Y. Savrasov, First-principles simulations of heavy Fermion cerium compounds based on the Kondo lattice, Phys. Rev. Lett. 103, 096403 (2009).
  * Haule and Kotliar [2009] K. Haule and G. Kotliar, Arrested Kondo effect and hidden order in URu2Si2, Nat. Phys. 5, 796 (2009).
  * Shinaoka  _et al._ [2015] H. Shinaoka, S. Hoshino, M. Troyer, and P. Werner, Phase diagram of pyrochlore iridates: All-in–all-out magnetic ordering and non-fermi-liquid properties, Phys. Rev. Lett. 115, 156401 (2015).
  * Mandal  _et al._ [2019] S. Mandal, K. Haule, K. M. Rabe, and D. Vanderbilt, Systematic beyond-DFT study of binary transition metal oxides, npj Comput. Mater. 5, 115 (2019).
  * Fujiwara  _et al._ [2022] H. Fujiwara, K. Terashima, J. Otsuki, N. Takemori, H. O. Jeschke, T. Wakita, Y. Yano, W. Hosoda, N. Kataoka, A. Teruya, M. Kakihana, M. Hedo, T. Nakama, Y. Ōnuki, K. Yaji, A. Harasawa, K. Kuroda, S. Shin, K. Horiba, H. Kumigashira, Y. Muraoka, and T. Yokoya, Anomalously large spin-dependent electron correlation in the nearly half-metallic ferromagnet ${\mathrm{CoS}}_{2}$, Phys. Rev. B 106, 085114 (2022).
  * Kuneš [2011] J. Kuneš, Efficient treatment of two-particle vertices in dynamical mean-field theory, Phys. Rev. B 83, 085102 (2011).
  * Yin  _et al._ [2014] Z. P. Yin, K. Haule, and G. Kotliar, Spin dynamics and orbital-antiphase pairing symmetry in iron-based superconductors, Nat. Phys. 10, 845 (2014).
  * Strand  _et al._ [2019] H. U. R. Strand, M. Zingl, N. Wentzell, O. Parcollet, and A. Georges, Magnetic response of ${\mathrm{Sr}}_{2}{\mathrm{RuO}}_{4}$: Quasi-local spin fluctuations due to Hund’s coupling, Phys. Rev. B 100, 125120 (2019).
  * Geffroy  _et al._ [2019] D. Geffroy, J. Kaufmann, A. Hariki, P. Gunacker, A. Hausoel, and J. Kuneš, Collective modes in excitonic magnets: Dynamical mean-field study, Phys. Rev. Lett. 122, 127601 (2019).
  * Shim  _et al._ [2007] J. H. Shim, K. Haule, and G. Kotliar, Modeling the Localized-to-Itinerant Electronic Transition in the Heavy Fermion System CeIrIn5, Science 318, 1615 (2007).
  * Otsuki  _et al._ [2009] J. Otsuki, H. Kusunose, and Y. Kuramoto, Evolution of a Large Fermi Surface in the Kondo Lattice, Phys. Rev. Lett. 102, 017202 (2009).
  * Otsuki [2015] J. Otsuki, Competing $d$-Wave and $p$-Wave Spin-Singlet Superconductivities in the Two-Dimensional Kondo Lattice, Phys. Rev. Lett. 115, 036404 (2015).
  * Kuramoto  _et al._ [2009] Y. Kuramoto, H. Kusunose, and A. Kiss, Multipole Orders and Fluctuations in Strongly Correlated Electron Systems, J. Phys. Soc. Jpn. 78, 072001 (2009).
  * Zirngiebl  _et al._ [1984] E. Zirngiebl, B. Hillebrands, S. Blumenröder, G. Güntherodt, M. Loewenhaupt, J. M. Carpenter, K. Winzer, and Z. Fisk, Crystal-field excitations in CeB6 studied by Raman and neutron spectroscopy, Phys. Rev. B 30, 4052 (1984).
  * Fujita  _et al._ [1980] T. Fujita, M. Suzuki, T. Komatsubara, S. Kunii, T. Kasuya, and T. Ohtsuka, Anomalous specific heat of CeB6, Solid State Commun. 35, 569 (1980).
  * Nakao  _et al._ [2001] H. Nakao, K.-i. Magishi, Y. Wakabayashi, Y. Murakami, K. Koyama, K. Hirota, Y. Endoh, and S. Kunii, Antiferro-Quadrupole Ordering of CeB6 Studied by Resonant X-Ray Scattering, J. Phys. Soc. Jpn. 70, 1857 (2001).
  * Tanaka  _et al._ [2004] Y. Tanaka, U. Staub, K. Katsumata, S. W. Lovesey, J. E. Lorenzo, Y. Narumi, V. Scagnoli, S. Shimomura, Y. Tabata, Y. Onuki, Y. Kuramoto, A. Kikkawa, T. Ishikawa, and H. Kitamura, Direct and quantitative determination of the orbital ordering in CeB6 by X-ray diffraction, Europhys. Lett. 68, 671 (2004).
  * Portnichenko  _et al._ [2020] P. Y. Portnichenko, A. Akbari, S. E. Nikitin, A. S. Cameron, A. V. Dukhnenko, V. B. Filipov, N. Y. Shitsevalova, P. Čermák, I. Radelytskyi, A. Schneidewind, J. Ollivier, A. Podlesnyak, Z. Huesges, J. Xu, A. Ivanov, Y. Sidis, S. Petit, J.-M. Mignot, P. Thalmeier, and D. S. Inosov, Field-angle-resolved magnetic excitations as a probe of hidden-order symmetry in CeB6, Phys. Rev. X 10, 021010 (2020).
  * Effantin  _et al._ [1985] J. Effantin, J. Rossat-Mignod, P. Burlet, H. Bartholin, S. Kunii, and T. Kasuya, Magnetic phase diagram of CeB6, J. Magn. Magn. Mater. 47-48, 145 (1985).
  * Tayama  _et al._ [1997] T. Tayama, T. Sakakibara, K. Tenya, H. Amitsuka, and S. Kunii, Magnetic Phase Diagram of CexLa1-xB6 Studied by Static Magnetization Measurement at Very Low Temperatures, J. Phys. Soc. Jpn. 66, 2268 (1997).
  * Kubo and Kuramoto [2004] K. Kubo and Y. Kuramoto, Octupole Ordering Model for the Phase IV of CexLa1-xB6, J. Phys. Soc. Jpn. 73, 216 (2004).
  * Mannix  _et al._ [2005] D. Mannix, Y. Tanaka, D. Carbone, N. Bernhoeft, and S. Kunii, Order parameter segregation in ${\mathrm{ce}}_{0.7}{\mathrm{la}}_{0.3}{\mathrm{b}}_{6}$: $4f$ octopole and $5d$ dipole magnetic order, Phys. Rev. Lett. 95, 117206 (2005).
  * Kusunose and Kuramoto [2005] H. Kusunose and Y. Kuramoto, Evidence for Octupole Order in Ce0.7La0.3B6 from Resonant X-ray Scattering, J. Phys. Soc. Jpn. 74, 3139 (2005).
  * Nakamura  _et al._ [1995] S. Nakamura, T. Goto, and S. Kunii, Magnetic phase diagrams of the dense kondo compounds CeB6 and Ce0.5La0.5B6, J. Phys. Soc. Jpn. 64, 3941 (1995).
  * J. Ohkawa [1985] F. J. Ohkawa, Orbital Antiferromagnetism in CeB6, J. Phys. Soc. Jpn. 54, 3909 (1985).
  * Ōnuki  _et al._ [1989] Y. Ōnuki, T. Komatsubara, P. H. P. Reinders, and M. Springford, Fermi Surface and Cyclotron Mass of CeB6, J. Phys. Soc. Jpn. 58, 3698 (1989).
  * Shiina  _et al._ [1997] R. Shiina, H. Shiba, and P. Thalmeier, Magnetic-Field Effects on Quadrupolar Ordering in a $\Gamma_{8}$-Quartet System CeB6, J. Phys. Soc. Jpn. 66, 1741 (1997).
  * Shiina  _et al._ [1998] R. Shiina, O. Sakai, H. Shiba, and P. Thalmeier, Interplay of Field-Induced Multipoles in CeB6, J. Phys. Soc. Jpn. 67, 941 (1998).
  * Shiba  _et al._ [1999] H. Shiba, O. Sakai, and R. Shiina, Nature of Ce-Ce Interaction in CeB6 and Its Consequences, J. Phys. Soc. Jpn. 68, 1988 (1999).
  * Sakurai and Kuramoto [2004] G. Sakurai and Y. Kuramoto, Multipolar Interactions in the Anderson Lattice with Orbital Degeneracy, J. Phys. Soc. Jpn. 73, 225 (2004).
  * Sakurai and Kuramoto [2005] G. Sakurai and Y. Kuramoto, Wavenumber Dependence of Multipolar Interactions in the Anderson Lattice, J. Phys. Soc. Jpn. 74, 975 (2005).
  * Yamada and Hanzawa [2019] T. Yamada and K. Hanzawa, Derivation of RKKY Interaction between Multipole Moments in CeB6 by the Effective Wannier Model Based on the Bandstructure Calculation, J. Phys. Soc. Jpn. 88, 084703 (2019).
  * Hanzawa and Yamada [2019] K. Hanzawa and T. Yamada, Origin of anisotropic RKKY interactions in CeB6, J. Phys. Soc. Jpn. 88, 124710 (2019).
  * Kubo and Hotta [2005] K. Kubo and T. Hotta, Microscopic theory of multipole ordering in NpO2, Phys. Rev. B 71, 140404 (2005).
  * Iizuka  _et al._ [2022] Y. Iizuka, T. Yamada, K. Hanzawa, and Y. Ōno, First-Principles Study of the RKKY Interaction and the Quadrupole Order in the Pr 1-2-20 Systems Pr$T_{2}$Al20 ($T$ = Ti, V), J. Phys. Soc. Jpn. 91, 074708 (2022).
  * Otsuki  _et al._ [2019] J. Otsuki, K. Yoshimi, H. Shinaoka, and Y. Nomura, Strong-coupling formula for momentum-dependent susceptibilities in dynamical mean-field theory, Phys. Rev. B 99, 165134 (2019).
  * Koepernik and Eschrig [1999] K. Koepernik and H. Eschrig, Full-potential nonorthogonal local-orbital minimum-basis band-structure scheme, Phys. Rev. B 59, 1743 (1999).
  * Tanaka  _et al._ [1997] K. Tanaka, Y. Kato, and Y. Onuki, $4f$-Electron Density Distribution in Crystals of CeB6 at 165 K and its Analysis Based on the Crystal Field Theory, Acta Crystallographica Section B 53, 143 (1997).
  * Eschrig and Koepernik [2009] H. Eschrig and K. Koepernik, Tight-binding models for the iron-based superconductors, Phys. Rev. B 80, 104503 (2009).
  * Koepernik  _et al._ [2021] K. Koepernik, O. Janson, Y. Sun, and J. van den Brink, Symmetry conserving maximally projected wannier functions (2021), arXiv:2111.09652 [cond-mat.mtrl-sci] .
  * Hotta [2004] T. Hotta, Spin and orbital structure of uranium compounds on the basis of a $j$-$j$ coupling scheme, Phys. Rev. B 70, 054405 (2004).
  * Hotta and Harima [2006] T. Hotta and H. Harima, Effective Crystalline Electric Field Potential in a $j$-$j$ Coupling Scheme, J. Phys. Soc. Jpn. 75, 124711 (2006).
  * Anisimov  _et al._ [1997] V. I. Anisimov, F. Aryasetiawan, and A. I. Lichtenstein, First-principles calculations of the electronic structure and spectra of strongly correlated systems: the LDA + U method, J. Phys.: Condens. Matter 9, 767 (1997).
  * Chiaia  _et al._ [1997] G. Chiaia, O. Tjernberg, L. Du‘o, S. De Rossi, P. Vavassori, I. Lindau, T. Takahashi, S. Kunii, T. Komatsubara, D. Cocco, S. Lizzit, and G. Paolucci, CeB6 and CeCu6 single crystals probed by resonant photoemission spectroscopy: A comparison between the two electronic structures, Phys. Rev. B 55, 9207 (1997).
  * Souma  _et al._ [2001] S. Souma, H. Kumigashira, T. Ito, T. Sato, T. Takahashi, and S. Kunii, Crystal-field splitting in CeB6 observed by ultrahigh-resolution photoemission spectroscopy, J. Electron Spectros. Relat. Phenomena 114-116, 729 (2001).
  * Lang  _et al._ [1981] J. K. Lang, Y. Baer, and P. A. Cox, Study of the 4f and valence band density of states in rare-earth metals. II. Experiment and results, J. Phys. F: Met. Phys. 11, 121 (1981).
  * Herbst  _et al._ [1976] J. F. Herbst, R. E. Watson, and J. W. Wilkins, $4f$ excitation energies in rare-earth metals: Relativistic calculations, Phys. Rev. B 13, 1439 (1976).
  * Herbst  _et al._ [1978] J. F. Herbst, R. E. Watson, and J. W. Wilkins, Relativistic calculations of $4f$ excitation energies in the rare-earth metals: Further results, Phys. Rev. B 17, 3089 (1978).
  * Shinaoka  _et al._ [2020a] H. Shinaoka, J. Otsuki, M. Kawamura, N. Takemori, and K. Yoshimi, DCore: Integrated DMFT software for correlated electrons, SciPost Phys. 10, 117 (2020a).
  * Parcollet  _et al._ [2015] O. Parcollet, M. Ferrero, T. Ayral, H. Hafermann, I. Krivenko, L. Messio, and P. Seth, TRIQS: A toolbox for research on interacting quantum systems, Comput. Phys. Commun. 196, 398 (2015).
  * Aichhorn  _et al._ [2016] M. Aichhorn, L. Pourovskii, P. Seth, V. Vildosola, M. Zingl, O. E. Peil, X. Deng, J. Mravlje, G. J. Kraberger, C. Martins, M. Ferrero, and O. Parcollet, TRIQS/DFTTools: A TRIQS application for ab initio calculations of correlated materials, Comput. Phys. Commun. 204, 200 (2016).
  * Antipov  _et al._ [2017] A. E. Antipov, I. Krivenko, and S. Iskakov, aeantipov/pomerol: 1.2 (2017).
  * Lebègue  _et al._ [2006a] S. Lebègue, A. Svane, M. I. Katsnelson, A. I. Lichtenstein, and O. Eriksson, Multiplet effects in the electronic structure of light rare-earth metals, Phys. Rev. B 74, 045114 (2006a).
  * Lebègue  _et al._ [2006b] S. Lebègue, A. Svane, M. I. Katsnelson, A. I. Lichtenstein, and O. Eriksson, Multiplet effects in the electronic structure of heavy rare-earth metals, J. Phys.: Condens. Matter 18, 6329 (2006b).
  * Locht  _et al._ [2016] I. L. M. Locht, Y. O. Kvashnin, D. C. M. Rodrigues, M. Pereiro, A. Bergman, L. Bergqvist, A. I. Lichtenstein, M. I. Katsnelson, A. Delin, A. B. Klautau, B. Johansson, I. Di Marco, and O. Eriksson, Standard model of the rare earths analyzed from the Hubbard I approximation, Phys. Rev. B 94, 085137 (2016).
  * Delange  _et al._ [2017] P. Delange, S. Biermann, T. Miyake, and L. Pourovskii, Crystal-field splittings in rare-earth-based hard magnets: An ab initio approach, Phys. Rev. B 96, 155132 (2017).
  * Harima  _et al._ [1988] H. Harima, O. Sakai, T. Kasuya, and A. Yanase, New interpretation of the de Haas-van Alphen signals of LaB6, Solid State Commun. 66, 603 (1988).
  * Koitzsch  _et al._ [2016] A. Koitzsch, N. Heming, M. Knupfer, B. Büchner, P. Y. Portnichenko, A. V. Dukhnenko, N. Y. Shitsevalova, V. B. Filipov, L. L. Lev, V. N. Strocov, J. Ollivier, and D. S. Inosov, Nesting-driven multipolar order in CeB6 from photoemission tomography, Nat. Commun. 7, 10876 (2016).
  * Jarrell [1992] M. Jarrell, Hubbard model in infinite dimensions: A quantum Monte Carlo study, Phys. Rev. Lett. 69, 168 (1992).
  * Kusunose [2008] H. Kusunose, Description of multipole in $f$-electron systems, J. Phys. Soc. Jpn. 77, 064710 (2008).
  * Gull  _et al._ [2011] E. Gull, A. J. Millis, A. I. Lichtenstein, A. N. Rubtsov, M. Troyer, and P. Werner, Continuous-time monte carlo methods for quantum impurity models, Rev. Mod. Phys. 83, 349 (2011).
  * Tagliavini  _et al._ [2018] A. Tagliavini, S. Hummel, N. Wentzell, S. Andergassen, A. Toschi, and G. Rohringer, Efficient Bethe-Salpeter equation treatment in dynamical mean-field theory, Phys. Rev. B 97, 235140 (2018).
  * Shinaoka  _et al._ [2020b] H. Shinaoka, D. Geffroy, M. Wallerberger, J. Otsuki, K. Yoshimi, E. Gull, and J. Kuneš, Sparse sampling and tensor network representation of two-particle Green’s functions, SciPost Phys. 8, 012 (2020b).
  * Wallerberger  _et al._ [2021] M. Wallerberger, H. Shinaoka, and A. Kauch, Solving the Bethe-Salpeter equation with exponential convergence, Phys. Rev. Research 3, 033168 (2021).
  * Kitatani  _et al._ [2022] M. Kitatani, R. Arita, T. Schäfer, and K. Held, Strongly correlated superconductivity with long-range spatial fluctuations, J. Phys. Mater. 5, 034005 (2022).