#  Material-based analysis of spin-orbital Mott insulators

Ryuta Iwazaki, Hiroshi Shinaoka, and Shintaro Hoshino  Department of Physics,
Saitama University, Shimo-Okubo, Saitama 338-8570, Japan

###### Abstract

We present a framework for analyzing Mott insulators using a material-based
tight-binding model. We start with a realistic multiorbital Hubbard model and
derive an effective model for the localized electrons through the second-order
perturbation theory with respect to intersite hopping. This effective model,
known as the Kugel-Khomskii model, is described by SU($N$) generators, where
$N$ is the number of localized states. We solve this model by the mean-field
theory that takes local correlations into account and reveal spin-orbital
ordered states. To include spatial correlations, we apply the classical Monte
Carlo based on the path-integral approach with SU($N$) coherent states, and
also derive the equation of motion for spin-orbital degrees of freedom. Our
approach is applicable to any Mott insulator with reasonable computational
cost. The $5d$-pyrochlore oxide is used here as demonstration.

Introduction.— Multiorbital systems with strongly correlated electrons have
been attracting attention due to their diverse physical phenomena, such as
electronic ordering and multiferroic behavior. It is crucial to uncover their
material-specific physical properties in order to make a serious comparison
with experimental results. In materials with weakly correlated electrons,
density functional theory (DFT)-based calculations have been successful in
describing their electronic properties. On the other hand, in the strongly
correlated regime, it is useful to construct a tight-binding model using
localized Wannier functions and subsequently employ a multiorbital Hubbard
model with local Coulomb repulsive interactions as a fundamental model.
Unfortunately, it is extremely difficult to perform the calculations in a
realistic setting due to the immense numerical cost. A theoretical framework
that is applicable to realistic strongly correlated electron systems is highly
desired, which will enable material prediction through, for example, high-
throughput screening [1].

In the present work, we focus on the Mott insulators where the electrons are
localized with strong local Coulomb interaction. Even in this case, the spin-
orbital degrees of freedom generate a number of interesting phenomena such as
magnetic orderings, multiferroic behaviors and spin liquids [2, 3, 4, 5, 6, 7,
8, 9, 10]. The low-energy effective model with localized electrons is known as
the Kugel-Khomskii model, in which both the spin and orbital degrees of
freedom are involved [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24,
25, 26, 27]. The realistic localized models have been discussed for the spin
model [28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39] and
$e_{g}$/$t_{2g}$-multiorbital systems [40, 41, 42, 43, 44, 45, 46]. The
DFT+DMFT approaches have also been employed for the analysis [40, 41, 28, 42,
47, 48, 39, 49]. In order to study arbitrary Mott insulator materials,
however, a more general framework is needed that can be applied at reasonable
computational cost to general multiorbital systems with spin-orbit
interactions and any number $N$ of localized states per atom.

In this paper, we propose a general framework to perform calculations for the
spin-orbital Mott insulators, which is not restricted to specific systems. We
develop a realistic Kugel-Khomskii model based on the tight-binding model
derived from the first-principles calculation and the local Coulomb
interaction with Slater-Condon parameters. The model contains $N^{2}-1$ spin-
orbital degrees of freedom and is described by SU($N$) generators. When
analyzing the model, while a fully quantum analysis is not feasible because of
a huge computational cost, we use the classical Monte Carlo with the SU($N$)
coherent state [50, 51, 52], in addition to the standard mean-field theory.
The SU($N$) coherent state has been used for the spin systems [53, 54, 55, 56,
57, 58, 59, 60, 61, 62], and here we apply it to the realistic Kugel-Khomskii
model. While the quantum mechanical inter-site correlations at very low
temperatures are not incorporated in our theory, our method captures the
characteristic physics at finite temperatures in a realistic setup for any
Mott insulators with reasonable numerical cost.

We will take the pyrochlore oxide Cd${}_{2}T_{2}$O7 as an example. This is
suitable as a prototype material for the demonstration of our framework due to
its complicated electronic structure: the four transition metal $T$ atoms in
unit cell (specified as sublattice indices A, B, C, D), large spin-orbit
interaction, and $t_{2g}$ three orbitals of $5d$ electrons with trigonal
symmetry at $T$ atom site [63, 64] (see Fig. 1). In addition, their non-
colinear magnetic structures are well studied both theoretically and
experimentally [65, 66]. Hence the applicability to this prototypical material
Cd${}_{2}T_{2}$O7 demonstrates the versatility of our method.

Realistic Kugel-Khomski model.— The realistic effective model for the
localized electrons are constructed based on the multiorbital Hubbard model
derived from the first principles calculation. Let us begin with the
Hamiltonian $\mathscr{H}=\mathscr{H}_{\mathrm{loc}}+\mathscr{H}_{t}$, where

$\displaystyle\mathscr{H}_{t}=\sum_{\langle
ij\rangle}\sum_{ab}t^{ab}_{ij}c_{ia}^{\dagger}c_{jb}+\mathrm{H.c.}$ (1)

describes the intersite hopping term. The operator $c_{ia}$ annihilates the
electron at the atom site $i$ with the spin($\sigma$)-orbital($\gamma$) index
$a=(\gamma,\sigma)$. The symbol $\langle ij\rangle$ indicates the summation
with respect to the pairs of atomic sites, and includes the terms other than
the nearest neighbor sites. The local part $\mathscr{H}_{\mathrm{loc}}$ is
further divided into three components as
$\mathscr{H}_{\mathrm{loc}}=\mathscr{H}_{U}+\mathscr{H}_{\mathrm{SOC}}+\mathscr{H}_{\mathrm{CEF}}$,
which are the Coulomb interaction, the spin-orbit coupling and the local
crystalline electric field, respectively. The Coulomb interaction is written
as

$\displaystyle\mathscr{H}_{U}=\sum_{i\gamma_{1}\gamma_{2}\gamma_{3}\gamma_{4}\sigma\sigma^{\prime}}U_{\gamma_{1}\gamma_{2}\gamma_{3}\gamma_{4}}c_{i\gamma_{1}\sigma}^{\dagger}c_{i\gamma_{2}\sigma^{\prime}}^{\dagger}c_{i\gamma_{4}\sigma^{\prime}}c_{i\gamma_{3}\sigma},$
(2)

which is parameterized by the Slater-Condon parameters as typically used in
LDA+$U$ or LDA+DMFT framework [67]. Specifically for the three orbital case as
in $t_{2g}$ orbital, the standard Slater-Kanamori form is employed:
$U_{\gamma\gamma\gamma\gamma}=U/2$,
$U_{\gamma\gamma^{\prime}\gamma\gamma^{\prime}}=U^{\prime}/2$,
$U_{\gamma\gamma^{\prime}\gamma^{\prime}\gamma}=U_{\gamma\gamma\gamma^{\prime}\gamma^{\prime}}=J/2$
for $\gamma\neq\gamma^{\prime}$ ($U^{\prime}=U-2J$) and the other terms are
zero.

In the following, we take the tight-binding model of $T=$ Os derived from
electronic-structure calculation [68]. Since the band structure [see Fig.
2(a)] is similar to the other materials with different filling such as $T=$ Re
[69, 70], we use the data of the $T=$ Os case also for the other electron
fillings.

Figure 1:  Single-site eigenenergy levels for Cd2$T_{2}$O7, which corresponds
to $\mathscr{H}_{t}=0$. We choose $U=4$ eV and $J/U=0.1$ which is comparable
to the previous study [71]. The vertical axis shows the energy measured from
the lowest energy at each $n$. The inset is the crystal structure of
Cd${}_{2}T_{2}$O7, where only $T$ and O atoms are shown for clarity [72].

We analyze the multiorbital Hubbard model in the strong coupling limit
($U\to\infty$), where the electrons are localized. First of all, we derive the
eigenenergies and eigenfunctions in the atomic model with only
$\mathscr{H}_{\rm loc}$, which is necessary for specifying the model Hilbert
space at low energies, i.e., the number $N$ of the localized states. Figure 1
shows the single-site eigenenergy diagram of $\mathscr{H}_{\mathrm{loc}}$ for
each number $n$ of electrons per $T$ atom. When we focus on the odd number of
the filling $n$, there are only doubly degenerated states corresponding to the
Kramers doublet.

In this paper, we choose $n=1$ for a demonstration of our scheme, which allows
us to choose the size of the model space as $N=2$, $4$, $6$ based on Fig. 1.
We call them SU(2), SU(4) and SU(6) models, respectively. The SU($N$) model
contains $N^{2}-1$ operators for each atom. The procedure for the simplest
$N=2$ case is summarized in Supplementary Matreial (SM) [68]. Although the
dimension of the model Hilbert space may be dependent on the lattice site, we
here take the same $N$ for all the sites.

Once the model space is specified, we treat the intersite Hamiltonian
$\mathscr{H}_{t}$ as a perturbation, to obtain the effective Hamiltonian which
gives a correct eigenenergies within the restricted Hilbert space [73, 74, 75,
76, 77]. While there are several choices of the form of the effective
Hamiltonian, the Hermitian Hamiltonian (des Cloizeaux type) is easier to be
handled [75, 76]. We focus on the two atoms which are connected by the hopping
matrix $\mathscr{H}_{t}$, and expand this two-site Hamiltonian up to second
order of $\mathscr{H}_{t}$ [68]. Thereby we obtain the matrix element of the
effective Hamiltonian whose size is $N^{2}\times N^{2}$. We can rewrite the
obtained effective Hamiltonian by complete local operators $\mathscr{O}_{i}$
at the site $i$. We employ the numerical calculation with matrix
multiplications for this procedure [78]. Collecting all the combinations of
the two-site Hamiltonians, we obtain the following realistic Kugel-Khomskii
model:

$\displaystyle\mathscr{H}_{\mathrm{eff}}[\mathscr{O}]=\sum_{\langle
ij\rangle}\sum_{\xi\xi^{\prime}}I_{ij}^{\xi\xi^{\prime}}\mathscr{O}_{i}^{\xi}\mathscr{O}_{j}^{\xi^{\prime}}-\sum_{i}\sum_{\xi}H_{i}^{\xi}\mathscr{O}_{i}^{\xi},$
(3)

where both the zeroth- and second-order contributions are involved in this
effective Hamiltonian. We have defined the local operators
$\mathscr{O}_{i}^{\xi}=\sum_{\alpha\beta}O_{\alpha\beta}^{\xi}|\alpha\rangle_{i}\,{}_{i}\langle\beta|$
($\alpha=1,\cdots,N$, $\xi=0,\cdots,N^{2}-1$), where $|\alpha\rangle_{i}$ is a
state vector in the model Hilbert space at site $i$. We use the matrices
$O^{\xi}_{\alpha\beta}$ with completeness and orthonormality (e.g. for single
orbital model, we take the SU(2) generators, which are the Pauli matrices)
[68]. We emphasize that this Hamiltonian is derived from the first-principles
calculation data, where the tunable parameters are only the local Coulomb
interaction parameters $U$ and $J$. In the actual calculation, the data of
$I_{ij}^{\xi\xi^{\prime}}$ is outputted with the data structure similar to the
original input of $t_{ij}^{ab}$.

Since it is in general difficult to interpret the physical meaning of the
local operators $\mathscr{O}_{i}^{\xi}$, it is desirable to transform them
into physical quantities defined in terms of the original electronic system.
Let us consider the local physical quantity $\mathscr{A}_{i}$. This can be
spin or orbital operator if we choose the form of
$\mathscr{A}_{i}=\frac{1}{2}\sum_{ab}A_{ab}c_{ia}^{\dagger}c_{ib}$ where the
matrix $A$ is composed of a direct product of the matrices in spin and orbital
spaces. By using the projection operator onto the model Hilbert space,
$\mathscr{P}=\prod_{i}\sum_{\alpha}|\alpha\rangle_{i}\,{}_{i}\langle\alpha|$,
we obtain

$\displaystyle\mathscr{P}\mathscr{A}_{i}\mathscr{P}$
$\displaystyle=\sum_{\xi}\mathscr{O}_{i}^{\xi}\sum_{\alpha\beta}{}_{i}\langle\alpha|\mathscr{A}_{i}|\beta\rangle_{i}O^{\xi}_{\beta\alpha}.$
(4)

We can get the matrix element
${}_{i}\langle\alpha|\mathscr{A}_{i}|\beta\rangle_{i}$ by analyzing
$\mathscr{H}_{\rm loc}$. Thus, once the expectation value of $\mathscr{O}_{i}$
is obtained by solving the model in Eq. (3), any local physical quantities can
be evaluated through this formula. It is notable that $\mathscr{A}_{i}$ can be
chosen as many-body quantities such as a double occupancy, which is not
usually considered for the conventional Kugel-Khomskii model.

The correlation functions are also useful quantities. When we consider the
linear response against a small fictitious field conjugate to
$\mathscr{O}_{i}^{\xi}$, the dynamical susceptibilities are given by

$\displaystyle\chi_{ij}^{\xi\xi^{\prime}}(\mathrm{i}\nu)=\int_{0}^{1/T}\mathrm{d}\tau\,\quantity[\langle\mathscr{O}_{i}^{\xi}(\tau)\mathscr{O}_{j}^{\xi^{\prime}}\rangle-\langle\mathscr{O}_{i}^{\xi}\rangle\langle\mathscr{O}_{j}^{\xi^{\prime}}\rangle]\mathrm{e}^{\mathrm{i}\nu\tau},$
(5)

where
$\mathscr{O}_{i}^{\xi}(\tau)=\mathrm{e}^{\tau\mathscr{H}}\mathscr{O}_{i}^{\xi}\mathrm{e}^{-\tau\mathscr{H}}$,
and $\tau$ is a Heisenberg picture with imaginary time, and $\nu=2\pi mT$
($m\in\mathbb{Z}$) is a bosonic Matsubara frequency. We have taken
$k_{\mathrm{B}}=1$. Using Eq. (4), the susceptibility can be transformed into
the physical susceptibilities defined in terms of the original electron
operators. The information of any spin-orbital excitation is encoded in Eq.
(5). For example, we can obtain the dispersion of the orbiton, which is a
quasiparticle describing the excitation of the orbital [13, 17, 19].

Mean-field theory.— Since the obtained localized model contains quantum
effects, it is still very hard to be solved. In the following, we introduce
several approximated methods to solve the realistic Kugel-Khomskii model given
in Eq. (3). The most fundamental approximation is the mean-field theory.
Defining the effective field $\tilde{H}_{i}^{\xi}=H_{i}^{\xi}-\sum_{j\neq
i,\xi^{\prime}}I_{ij}^{\xi\xi^{\prime}}M_{j}^{\xi^{\prime}}$, the mean-field
Hamiltonian is written as

$\displaystyle\mathscr{H}_{\mathrm{MF}}=-\sum_{i}\sum_{\xi}\tilde{H}_{i}^{\xi}\mathscr{O}_{i}^{\xi}-\sum_{\langle
ij\rangle}\sum_{\xi\xi^{\prime}}I_{ij}^{\xi\xi^{\prime}}M_{i}^{\xi}M_{j}^{\xi^{\prime}}.$
(6)

We have defined
$M_{i}^{\xi}=\langle\mathscr{O}_{i}^{\xi}\rangle_{\mathrm{MF}}$ where the
expectation value is taken as
$\langle\cdots\rangle_{\mathrm{MF}}=\Tr(\cdots\mathrm{e}^{-\mathscr{H}_{\mathrm{MF}}/T})/\Tr\mathrm{e}^{-\mathscr{H}_{\mathrm{MF}}/T}$.
We also evaluate the dynamical susceptibilities with the random phase
approximation as

$\displaystyle\hat{\chi}(\bm{q},\omega)=\hat{\chi}_{0}(\omega)\quantity[\hat{1}+\hat{I}(\bm{q})\hat{\chi}_{0}(\omega)]^{-1},$
(7)

where the hat ($\,\hat{}\,$) symbol represents the matrix with respect to the
index $\xi$, and $\hat{1}$ is the identity matrix. We have defined the local
susceptibility by
$\hat{\chi}_{0}(\omega)=\hat{\chi}_{ii}(\omega+\mathrm{i}0^{+})$ which is
evaluated by the local mean-field Hamiltonian.

Figure 2:  (a) Electronic energy band for Cd${}_{2}T_{2}$O7. The vertical axis
is measured from the bottom of the bands. The horizontal dashed lines express
the chemical potential for each $n$. (b) Bosonic energy spectra
$\imaginary\chi(\bm{q},\omega)/\omega$ for the SU(6) model at $T=10^{-3}$ eV.
The left panel shows the dispersion for the spin, while the orbital excitation
spectra is shown in the right panel.

First of all, we show in Fig. 2(b) the spin-orbital excitation spectra of the
realistic Kugel-Khomskii model, which is contrasted against the fermionic
excitation of the original tight-binding electrons in (a). We take the SU(6)
model at $n=1$ and $T=10^{-3}$ eV. The left panel of (b) is the spectra for
the spin, which corresponds to the dispersion of the magnon. The gapped
excitation reflects the presence of the spin-orbit coupling. The right panel
is the spectra for the non-magnetic orbital (quadrupole) moment (see Ref. [68]
for the definition of the orbital moment). This orbital excitation is unique
to the SU(6) model, although the magnon dispersion is captured already in the
SU(2) model.

Figure 3:  Temperature dependence of (a) the order parameters and (b) the
$\bm{q}=\bm{0}$ component of the static diagonal susceptibilities at A-site
for the SU(6) model obtained by the mean-field analysis. The blue, orange,
green, red and purple lines show the spin, magnetic orbital, electric orbital,
electric dipole and magnetic octupole, respectively. (c) Sketches for the spin
($S$) and the magnetic orbital moment ($L$) of the SU(6) model at $T=10^{-3}$
eV. (d) Temperature dependence of the entropy (upper panel) and the specific
heat (lower panel) with the circles, crosses and triangles for the SU(2),
SU(4) and SU(6) models, respectively. For clarity, specific heat is vertically
shifted for the SU(4,6) models.

We show the temperature dependence of the order parameters at A-sublattice in
Fig. 3(a) for the SU(6) model. The symbols $S,L,Q,G$ and $T$ are the spin,
magnetic orbital, electric orbital (quadrupole), electric dipole and magnetic
octupole moments, respectively [68]. At low temperatures with $T\lesssim
10^{-2}$ eV, the magnetic ordering occurs, whose order parameters are
described by $S,L$ and $T$. We also show the $\bm{q}=\bm{0}$ component of the
diagonal susceptibilities at A-sublattice in Fig. 3(b), where the magnetic
susceptibilities ($S$, $L$, $T$) diverge. The magnetic structures at
$T=10^{-3}$ eV are shown at Fig. 3(c), which displays the all-in-all-out
(AIAO) structure and the antiparallel alignment of $S$ and $L$ moment. The
AIAO-type magnetic ordering in $5d$ pyrochlore oxides have been suggested both
theoretically and experimentally [79, 65, 66, 80, 81, 71].

The temperature dependence of the thermodynamic quantities per site are shown
in Fig. 3(d), with which we compare the results of the SU(2,4,6) models. All
of the models have an anomaly in the specific heat (lower panel) around
$T_{\mathrm{c}}\simeq 10^{-2}$ eV, which signals a second-order phase
transition. The SU(6) model has a smaller magnetic transition temperature
compared to the SU(2,4) cases. The single site entropy (upper panel) has a
$\ln 2$ plateau for the SU(4) and the SU(6) model just above $T_{\mathrm{c}}$,
and it deviates from $\ln 2$ reflecting the additional degrees of freedom at
higher $T$. The specific heat above $T_{\mathrm{c}}$ shows Schottky peaks
originating from the local energy-level splitting.

Classical model.—We can also solve the model by applying the classical
approximation to Eq. (3). In this method, we can examine the effect of the
non-local correlation. We employ the path-integral formalism using a coherent
state [52, 55], with which we derive both the classical partition function and
equations of motion. The coherent state is defined for each site $i$ by

$\displaystyle|\Omega_{i}\rangle=\sum_{\alpha=1}^{N}c_{\alpha}(\Omega_{i})|\alpha\rangle_{i},$
(8)

where $|\alpha\rangle_{i}$ is a quantum state basis. $\Omega_{i}$ is a set of
local continuous variables:
$\Omega_{i}=\\{\xi_{1i},\cdots,\xi_{N-1,i},\varphi_{1i},\cdots,\varphi_{N-1,i}\\}$,
each of which is written as $\Omega_{pi}$ ($p=1,\cdots,2(N-1)$) [68]. Here
$\xi_{1i,\cdots}\in[0,\pi/2]$ and $\varphi_{1i,\cdots}\in[0,2\pi)$
respectively correspond to the generalized versions of polar angle and
azimuthal angle of the spin in the SU(2) model.

The partition function is written as
$Z=\int\mathscr{D}\bm{\Omega}\,\mathrm{e}^{-\mathscr{S}}$, where the action is
[53]

$\displaystyle\mathscr{S}$
$\displaystyle=\int\mathrm{d}\tau\,\quantity(\langle\bm{\Omega}|\partial_{\tau}|\bm{\Omega}\rangle+\langle\bm{\Omega}|\mathscr{H}_{\rm
eff}|\bm{\Omega}\rangle).$ (9)

We have defined $|\bm{\Omega}\rangle=\prod_{i}|\Omega_{i}\rangle$ at an
imaginary time $\tau$. The quantum-mechanical operator $\mathscr{O}_{i}$ is
now replaced by the classical variable:
$\mathcal{O}^{\xi}(\Omega_{i})=\langle\bm{\Omega}|\mathscr{O}_{i}^{\xi}|\bm{\Omega}\rangle$.
Based on these expressions, the classical model can be rigorously derived by
using the coherent state path integral method that omits the Berry phase term,
as in the spin model [82]. We can also show that the classical free energy is
always larger than the quantum one [83], and it is ensured that the lowest-
free-energy state in the classical model is energetically closest to the
genuine quantum state.

The model can be numerically simulated by using the classical Monte Carlo
method. We use the local Metropolis update and the replica exchange method
which allow us to simulate the systems with various temperatures efficiently
[84]. In addition, we also apply the over-relaxation update [85] for the more
efficient simulation. The over-relaxation update in the present case consists
of microcanonical moves that does not alter the energy. For the SU(2) case,
the local spin vector is rotated around the local effective field by the angle
$\pi$ [86, 87]. However, this cannot be directly extended to SU($N$) case, and
the consideration based on the coherent state is needed.

To perform the over-relaxation update for the SU($N$) case, let us focus on
the one lattice site $i$, and then its effective local Hamiltonian is written
as $\mathcal{H}_{{\rm
loc},i}=-\sum_{\xi}\tilde{H}_{i}^{\xi}\mathcal{O}^{\xi}(\Omega_{i})$ where the
effect of the surrounding sites is included in
$\tilde{H}_{i}^{\xi}=H_{i}^{\xi}-\sum_{j\neq
i,\xi^{\prime}}I_{ij}^{\xi\xi^{\prime}}\mathcal{O}(\Omega_{j})$, which is not
dependent on $\Omega_{i}$. We can cast it into the coherent state
representation as

$\displaystyle\mathcal{H}_{{\rm loc},i}$
$\displaystyle=\sum_{\alpha\beta}h_{\alpha\beta}(i)c_{\alpha}^{*}(\Omega_{i})c_{\beta}(\Omega_{i})=\sum_{\gamma}\Lambda_{\gamma}(i)|d_{\gamma}(\Omega_{i})|^{2},$
(10)

where the diagonalization is performed in the right-most side by the unitary
matrix $V$: $d_{\gamma}=\sum_{\alpha}V^{\dagger}_{\gamma\alpha}c_{\alpha}$. It
is apparent at this point that the energy does not change by the phase
transformation $d_{\gamma}\to
d_{\gamma}\mathrm{e}^{\mathrm{i}\theta_{\gamma}}$ for any $\theta_{\gamma}$,
with which the coherent state is transformed as
$\Omega_{i}\to\Omega_{i}^{\prime}$. The parameter $\theta_{\gamma}$ is
determined to minimize the norm of the inner product
$\langle\Omega_{i}|\Omega_{i}^{\prime}\rangle$ (see Ref. [68] for more
details). This update makes it efficient to sample different configurations.
We note that the above procedure involving coherent state reproduces the over-
relaxation update usually used for the SU(2) case.

Figure 4:  (a) Temperature dependence of the specific heat per site for the
SU(2,4,6) models indicated by circles, crosses and triangles, respectively,
which are obtained by classical Monte Carlo method. The black dotted lines for
the SU(2,4) models are the results with $N_{\rm site}=256$ ($=4\times 4^{3}$).
(b) Temperature dependence of the $\bm{q}=\bm{0}$ component of the diagonal
susceptibilities at A-site for the SU(6) model.

We show the numerical result of the classical Monte Carlo in Fig. 4. The
calculation is performed for a finite-sized lattice with $N_{\rm site}=108$
($=4\times 3^{3}$) atoms where the lattice is created using primitive
translation vectors. The temperature dependence of the specific heat is shown
in Fig. 4(a) for the SU(2,4,6) models. At low temperatures, the specific heat
takes $2(N-1)\times\frac{1}{2}$ for the SU($N$) model, which satisfies the
equipartition theorem. Compared to the corresponding results of the mean-field
calculation Fig. 3(d), every model has the suppressed transition temperatures
down to $T_{\mathrm{c}}\sim 10^{-3}$ eV because of the incorporation of
spatial fluctuations. Figure 4(b) shows the $\bm{q}=\bm{0}$ component of the
static susceptibilities for the SU(6) model. For $T\lesssim 10^{-3}$ eV, each
component of the magnetic moments takes the huge values, showing the feature
of second-order phase transition. The electric ($Q,G$) susceptibilities are
characteristic for the SU(6) model and is absent in SU(2) cases.

Classical equation of motion.— Using the path-integral approach, our framework
can further address the thermodynamic non-equilibrium state. The equation of
motion itself has already been derived by Zhang-Batista [55]. The derivation
is based on the Heisenberg equation of motion of $\mathscr{O}_{i}^{\xi}$ which
gives $N^{2}-1$ equations. In terms of the parameters of the coherent states,
on the other hand, we only need $2(N-1)$ equations. Hence some of those
equations should be redundant. Here, we derive the $2(N-1)$ equations based on
the principle of the least action of Eq. (9) [82, 88]. The resultant equation
of motion for the local variable is given by

$\displaystyle\sum_{q}\mathcal{B}_{pq}(i)\frac{\partial\Omega_{qi}}{\partial\tau}=-\frac{\partial\mathcal{H}}{\partial\Omega_{pi}},$
(11)

where $\mathcal{H}=\langle\bm{\Omega}|\mathscr{H}_{\rm
eff}|\bm{\Omega}\rangle$ and the Berry curvature matrix is defined by

$\displaystyle\mathcal{B}_{pq}(i)$
$\displaystyle=\sum_{\alpha}\left(\frac{\partial
c_{\alpha}^{*}(\Omega_{i})}{\partial\Omega_{pi}}\frac{\partial
c_{\alpha}(\Omega_{i})}{\partial\Omega_{qi}}-\frac{\partial
c_{\alpha}^{*}(\Omega_{i})}{\partial\Omega_{qi}}\frac{\partial
c_{\alpha}(\Omega_{i})}{\partial\Omega_{pi}}\right),$ (12)

with $p,q=1,\cdots,2(N-1)$. Changing the time variable as
$\tau\to\mathrm{i}t$, we obtain the real-time equation of motion. Since the
analytic form of the Berry curvature matrix is obtained once the specific
coherent state is given in Eq. (8), the even-dimension antisymmetric matrix
$\mathcal{B}$ in Eq. (11) is easily inverted numerically. Thus the explicit
equation of motion is obtained for the $2(N-1)$ classical variables, and will
be used for a non-equilibrium dynamics in a realistic setup. The relation to
the equations in Ref. [55] is not apparent but can be deduced from the
equation

$\displaystyle-\frac{\partial\mathcal{O}_{i}^{\xi}}{\partial\tau}=\sum_{pq}\mathcal{B}^{-1}_{pq}(i)\frac{\partial\mathcal{O}_{i}^{\xi}}{\partial\Omega_{pi}}\frac{\partial\mathcal{H}}{\partial\Omega_{qi}},$
(13)

which derives from Eq. (11). The right-hand side is reminiscent of the
commutator $[\mathscr{O}_{i}^{\xi},\mathscr{H}_{\rm eff}]$.

Summary and outlook.— We have proposed the numerical calculation method for
generic spin-orbital Mott insulators, and applied it to $5d$-pyrochlore oxides
as a demonstration. A detailed comparison between simulation results and
experiments will provide us a deeper understanding of the Mott insulators,
which leads to a design of functional materials.

## Acknowledgement

The authors thank R. Pohle for fruitful discussions. This work was supported
by KAKENHI Grants No. 19H01842, No. 21K03459 and No. JP22J10620.

## References

  * Jain  _et al._ [2016] A. Jain, Y. Shin, and K. A. Persson, Nat. Rev. Mater. 1, 15004 (2016).
  * Imada  _et al._ [1998] M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys. 70, 1039 (1998).
  * Tokura and Nagaosa [2000] Y. Tokura and N. Nagaosa, Science 288, 462 (2000), doi: 10.1126/science.288.5465.462.
  * Kim  _et al._ [2008] B. J. Kim, H. Jin, S. J. Moon, J.-Y. Kim, B.-G. Park, C. S. Leem, J. Yu, T. W. Noh, C. Kim, S.-J. Oh, J.-H. Park, V. Durairaj, G. Cao, and E. Rotenberg, Phys. Rev. Lett. 101, 076402 (2008).
  * Jackeli and Khaliullin [2009] G. Jackeli and G. Khaliullin, Phys. Rev. Lett. 102, 017205 (2009).
  * Kim  _et al._ [2009] B. J. Kim, H. Ohsumi, T. Komesu, S. Sakai, T. Morita, H. Takagi, and T. Arima, Science 323, 1329 (2009).
  * Witczak-Krempa  _et al._ [2014] W. Witczak-Krempa, G. Chen, Y. B. Kim, and L. Balents, Annu. Rev. Condens. Matter Phys. 5, 57 (2014), doi: 10.1146/annurev-conmatphys-020911-125138.
  * Kitagawa  _et al._ [2018] K. Kitagawa, T. Takayama, Y. Matsumoto, A. Kato, R. Takano, Y. Kishimoto, S. Bette, R. Dinnebier, G. Jackeli, and H. Takagi, Nature 554, 341 (2018).
  * Kasahara  _et al._ [2018] Y. Kasahara, T. Ohnishi, Y. Mizukami, O. Tanaka, S. Ma, K. Sugii, N. Kurita, H. Tanaka, J. Nasu, Y. Motome, T. Shibauchi, and Y. Matsuda, Nature 559, 227 (2018).
  * Tang  _et al._ [2022] N. Tang, Y. Gritsenko, K. Kimura, S. Bhattacharjee, A. Sakai, M. Fu, H. Takeda, H. Man, K. Sugawara, Y. Matsumoto, Y. Shimura, J. Wen, C. Broholm, H. Sawa, M. Takigawa, T. Sakakibara, S. Zherlitsyn, J. Wosnitza, R. Moessner, and S. Nakatsuji, Nat. Phys. 10.1038/s41567-022-01816-4 (2022).
  * Kugel and Khomskii [1972] K. I. Kugel and D. I. Khomskii, Pis’ma Zh. Eksp. Teor. Fiz. 15, 629 (1972), [Sov. Phys. JETP Lett. 15, 446 (1972)].
  * Kugel and Khomskii [1973] K. I. Kugel and D. I. Khomskii, Zh. Eksp. Teor. Fiz. 64, 1429 (1973), [Sov. Phys. JETP 37, 725 (1973)].
  * Cyrot and Lyon-Caen [1975] M. Cyrot and C. Lyon-Caen, Le J. de Phys. 36, 253 (1975).
  * Khaliullin and Oudovenko [1997] G. Khaliullin and V. Oudovenko, Phys. Rev. B 56, R14243 (1997).
  * Ishihara  _et al._ [1997] S. Ishihara, J. Inoue, and S. Maekawa, Phys. Rev. B 55, 8280 (1997).
  * Feiner and Oleś [1999] L. Feiner and A. Oleś, Phys. Rev. B 59, 3295 (1999).
  * Ishihara and Maekawa [2000] S. Ishihara and S. Maekawa, Phys. Rev. B 62, 2338 (2000).
  * Harris  _et al._ [2004] A. B. Harris, A. Aharony, O. Entin-Wohlman, I. Y. Korenblit, and T. Yildirim, Phys. Rev. B 69, 094409 (2004).
  * Ishihara [2004] S. Ishihara, Phys. Rev. B 69, 075118 (2004).
  * Normand and Oleś [2008] B. Normand and A. M. Oleś, Phys. Rev. B 78, 094427 (2008).
  * Nasu and Ishihara [2013] J. Nasu and S. Ishihara, Phys. Rev. B 88, 205110 (2013).
  * Koga  _et al._ [2018] A. Koga, S. Nakauchi, and J. Nasu, Phys. Rev. B 97, 094427 (2018).
  * Otsuki  _et al._ [2019] J. Otsuki, K. Yoshimi, H. Shinaoka, and Y. Nomura, Phys. Rev. B 99, 165134 (2019).
  * Bieniasz  _et al._ [2019] K. Bieniasz, M. Berciu, and A. M. Oleś, Phys. Rev. B 100, 125109 (2019).
  * Nasu and Naka [2021] J. Nasu and M. Naka, Phys. Rev. B 103, L121104 (2021).
  * Khaliullin  _et al._ [2021] G. Khaliullin, D. Churchill, P. P. Stavropoulos, and H.-Y. Kee, Phys. Rev. Research 3, 033163 (2021).
  * Khomskii [2022] D. I. Khomskii, ECS J. Solid State Sci. Technol. 11, 054004 (2022).
  * Zhang  _et al._ [2012] G. Zhang, E. Gorelov, E. Koch, and E. Pavarini, Phys. Rev. B 86, 184413 (2012).
  * Chiesa  _et al._ [2013] A. Chiesa, S. Carretta, P. Santini, G. Amoretti, and E. Pavarini, Phys. Rev. Lett. 110, 157204 (2013).
  * Yamaji  _et al._ [2014] Y. Yamaji, Y. Nomura, M. Kurita, R. Arita, and M. Imada, Phys. Rev. Lett. 113, 107201 (2014).
  * Rau  _et al._ [2014] J. G. Rau, E. K.-H. Lee, and H.-Y. Kee, Phys. Rev. Lett. 112, 077204 (2014).
  * Winter  _et al._ [2016] S. M. Winter, Y. Li, H. O. Jeschke, and R. Valentí, Phys. Rev. B 93, 214431 (2016).
  * Winter  _et al._ [2017] S. M. Winter, K. Riedl, and R. Valentí, Phys. Rev. B 95, 060404 (2017).
  * Kurzydłowski and Grochala [2017] D. Kurzydłowski and W. Grochala, Phys. Rev. B 96, 155140 (2017).
  * Chiesa  _et al._ [2019] A. Chiesa, E. Macaluso, P. Santini, S. Carretta, and E. Pavarini, Phys. Rev. B 99, 235145 (2019).
  * Huang  _et al._ [2020] Z. Huang, D. Liu, A. Mansikkamäki, V. Vieru, N. Iwahara, and L. F. Chibotaru, Phys. Rev. Research 2, 033430 (2020).
  * Kaib  _et al._ [2021] D. A. S. Kaib, S. Biswas, K. Riedl, S. M. Winter, and R. Valentí, Phys. Rev. B 103, L140402 (2021).
  * Churchill and Kee [2022] D. Churchill and H.-Y. Kee, Phys. Rev. B 105, 014438 (2022).
  * Mosca  _et al._ [2022] D. F. Mosca, L. V. Pourovskii, and C. Franchini, Phys. Rev. B 106, 035127 (2022).
  * Pavarini  _et al._ [2008] E. Pavarini, E. Koch, and A. I. Lichtenstein, Phys. Rev. Lett. 101, 266405 (2008).
  * Pavarini and Koch [2010] E. Pavarini and E. Koch, Phys. Rev. Lett. 104, 086402 (2010).
  * Autieri  _et al._ [2014] C. Autieri, E. Koch, and E. Pavarini, Phys. Rev. B 89, 155109 (2014).
  * Snamina and Oleś [2016] M. Snamina and A. M. Oleś, Phys. Rev. B 94, 214426 (2016).
  * Jeanneau  _et al._ [2017] J. Jeanneau, P. Toulemonde, G. Remenyi, A. Sulpice, C. Colin, V. Nassif, E. Suard, E. S. Colera, G. R. Castro, F. Gay, C. Urdaniz, R. Weht, C. Fevrier, A. Ralko, C. Lacroix, A. A. Aligia, and M. Núñez-Regueiro, Phys. Rev. Lett. 118, 207207 (2017).
  * Aligia and Helman [2019] A. A. Aligia and C. Helman, Phys. Rev. B 99, 195150 (2019).
  * Zhang  _et al._ [2022a] X.-J. Zhang, E. Koch, and E. Pavarini, Phys. Rev. B 105, 115104 (2022a).
  * Pourovskii  _et al._ [2021] L. V. Pourovskii, D. F. Mosca, and C. Franchini, Phys. Rev. Lett. 127, 237201 (2021).
  * Mosca  _et al._ [2021] D. F. Mosca, L. V. Pourovskii, B. H. Kim, P. Liu, S. Sanna, F. Boscherini, S. Khmelevskyi, and C. Franchini, Phys. Rev. B 103, 104401 (2021).
  * Pourovskii [2022] L. V. Pourovskii, arXiv:2212.03358 (2022).
  * Perelomov [1972] A. M. Perelomov, Commun. Math. Phys. 26, 222 (1972).
  * Gnutzmann and Kus [1998] S. Gnutzmann and M. Kus, J. Phys. A: Math. Gen. 31, 9871 (1998).
  * Nemoto [2000] K. Nemoto, J. Phys. A: Math. Gen. 33, 3493 (2000).
  * Read and Sachdev [1989] N. Read and S. Sachdev, Nucl. Phys. B 316, 609 (1989).
  * Stoudenmire  _et al._ [2009] E. M. Stoudenmire, S. Trebst, and L. Balents, Phys. Rev. B 79, 214436 (2009).
  * Zhang and Batista [2021] H. Zhang and C. D. Batista, Phys. Rev. B 104, 104409 (2021).
  * Remund  _et al._ [2022] K. Remund, R. Pohle, Y. Akagi, J. Romhányi, and N. Shannon, Phys. Rev. Research 4, 033106 (2022).
  * Dahlbom  _et al._ [2022a] D. Dahlbom, H. Zhang, C. Miles, X. Bai, C. D. Batista, and K. Barros, Phys. Rev. B 106, 054423 (2022a).
  * Seifert and Savary [2022] U. F. P. Seifert and L. Savary, Phys. Rev. B 106, 195147 (2022).
  * Dahlbom  _et al._ [2022b] D. Dahlbom, C. Miles, H. Zhang, C. D. Batista, and K. Barros, Phys. Rev. B 106, 235154 (2022b).
  * Zhang  _et al._ [2022b] H. Zhang, Z. Wang, D. Dahlbom, K. Barros, and C. D. Batista, arXiv:2203.15248 (2022b).
  * Do  _et al._ [2022] S.-H. Do, H. Zhang, D. A. Dahlbom, T. J. Williams, V. O. Garlea, T. Hong, T.-H. Jang, S.-W. Cheong, J.-H. Park, K. Barros, C. D. Batista, and A. D. Christianson, arXiv:2205.11770 (2022).
  * Pohle  _et al._ [2022] R. Pohle, N. Shannon, and Y. Motome, arXiv:2212.10040 (2022).
  * Gardner  _et al._ [2010] J. S. Gardner, M. J. P. Gingras, and J. E. Greedan, Rev. Mod. Phys. 82, 53 (2010).
  * Shinaoka  _et al._ [2019] H. Shinaoka, Y. Motome, T. Miyake, S. Ishibashi, and P. Werner, J. Phys.: Condens. Matter 31, 323001 (2019).
  * Shinaoka  _et al._ [2012] H. Shinaoka, T. Miyake, and S. Ishibashi, Phys. Rev. Lett. 108, 247204 (2012).
  * Yamaura  _et al._ [2012] J. Yamaura, K. Ohgushi, H. Ohsumi, T. Hasegawa, I. Yamauchi, K. Sugimoto, S. Takeshita, A. Tokuda, M. Takata, M. Udagawa, M. Takigawa, H. Harima, T. Arima, and Z. Hiroi, Phys. Rev. Lett. 108, 247205 (2012).
  * Kotliar  _et al._ [2006] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Rev. Mod. Phys. 78, 865 (2006).
  * [68] See Supplementary Material.
  * Singh  _et al._ [2002] D. J. Singh, P. Blaha, K. Schwarz, and J. O. Sofo, Phys. Rev. B 65, 155109 (2002).
  * Harima [2002] H. Harima, J. Phys. Chem. Solids 63, 1035 (2002).
  * Shinaoka  _et al._ [2015] H. Shinaoka, S. Hoshino, M. Troyer, and P. Werner, Phys. Rev. Lett. 115, 156401 (2015).
  * Momma and Izumi [2011] K. Momma and F. Izumi, J. Appl. Crystallogr. 44, 1272 (2011).
  * Ôkubo [1954] S. Ôkubo, Prog. Theor. Phys. 12, 603 (1954).
  * Bloch [1958] C. Bloch, Nucl. Phys. 6, 329 (1958).
  * des Cloizeaux [1960] J. des Cloizeaux, Nucl. Phys. 20, 321 (1960).
  * Durand [1983] P. Durand, Phys. Rev. A 28, 3184 (1983).
  * Kuramoto [2020] Y. Kuramoto,  _Quantum Many-Body Physics_, Vol. 934 (Springer Japan, 2020).
  * Iwazaki and Hoshino [2021] R. Iwazaki and S. Hoshino, Phys. Rev. B 103, 235145 (2021).
  * Tomiyasu  _et al._ [2012] K. Tomiyasu, K. Matsuhira, K. Iwasa, M. Watahiki, S. Takagi, M. Wakeshima, Y. Hinatsu, M. Yokoyama, K. Ohoyama, and K. Yamada, J. Phys. Soc. Jpn. 81, 034709 (2012), doi: 10.1143/JPSJ.81.034709.
  * Sagayama  _et al._ [2013] H. Sagayama, D. Uematsu, T. Arima, K. Sugimoto, J. J. Ishikawa, E. O’Farrell, and S. Nakatsuji, Phys. Rev. B 87, 100403 (2013).
  * Disseler [2014] S. M. Disseler, Phys. Rev. B 89, 140413 (2014).
  * Auerbach [1994] A. Auerbach,  _Interacting Electrons and Quantum Magnetism_ (Springer New York, 1994).
  * Lieb [1973] E. H. Lieb, Commun. Math. Phys. 31, 327 (1973).
  * Hukushima and Nemoto [1996] K. Hukushima and K. Nemoto, J. Phys. Soc. Jpn. 65, 1604 (1996).
  * Creutz [1987] M. Creutz, Phys. Rev. D 36, 515 (1987).
  * Landau and Binder [2021] D. Landau and K. Binder,  _A guide to Monte Carlo simulations in statistical physics_ (Cambridge university press, 2021).
  * Alonso  _et al._ [1996] J. L. Alonso, A. Tarancón, H. G. Ballesteros, L. A. Fernández, V. Martín-Mayor, and A. M. Sudupe, Phys. Rev. B 53, 2537 (1996).
  * Nagaosa [1999] N. Nagaosa,  _Quantum Field Theory in Condensed Matter Physics_ (Springer Berlin Heidelberg, 1999).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, J. Phys.: Condens. Matter 21, 395502 (2009).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, J. Phys.: Condens. Matter 29, 465901 (2017).
  * Mostofi  _et al._ [2014] A. A. Mostofi, J. R. Yates, G. Pizzi, Y.-S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Comput. Phys. Commun. 185, 2309 (2014).
  * Dal Corso [2014] A. Dal Corso, Comput. Mater. Sci. 95, 337 (2014).
  * Blöchl [1994] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Mandrus  _et al._ [2001] D. Mandrus, J. R. Thompson, R. Gaal, L. Forro, J. C. Bryan, B. C. Chakoumakos, L. M. Woods, B. C. Sales, R. S. Fishman, and V. Keppens, Phys. Rev. B 63, 195104 (2001).
  * Klein [1974] D. J. Klein, J. Chem. Phys. 61, 786 (1974).
  * Georgi [2000] H. Georgi,  _Lie Algebras in Particle Physics_ , 2nd ed. (Perseus Book, 2000).

SUPPLEMENTARY MATERIAL FOR
“Material-based analysis of spin-orbital Mott insulators”
R. Iwazaki, H. Shinaoka, and S. Hoshino
(Dated: )

## SM 1: Details of the first-principles calculation

For constructing the tight-binding Hamiltonian, we used Quantum ESPRESSO [89,
90] and wannier90 [91]. In the band calculations using Quantum ESPRESSO, we
used pseudopotentials from pslibrary 1.0.0 [92] and a kinetic energy cutoff of
75 Ry for the PAW method [93]. The band calculations were done with the
experimental lattice structure at 180 K: $a=10.1598$ Å and
$x(\mathrm{O}_{1})=0.319$ [94]. We constructed maximally localized Wannier
functions using wannier90 for the $t_{\mathrm{2g}}$ manifold.

## SM 2: Effective Hamiltonian

We write the Hamiltonian as $H=H_{0}+V$ where $V$ is treated by the
perturbation theory. The Hermitian effective Hamiltonian is given by [75, 95,
76]

$\displaystyle\mathscr{H}_{\rm eff}$
$\displaystyle=(\Omega^{\dagger}\Omega)^{-1/2}\Omega^{\dagger}H\Omega(\Omega^{\dagger}\Omega)^{-1/2},$
(S1)

where $\Omega$ is the wave operator determined by the operator equation

$\displaystyle[H_{0},\Omega]=-V\Omega+\Omega V\Omega.$ (S2)

We introduce the projection operators $P_{0}$ onto the model Hilbert space and
also $Q_{0}=1-P_{0}$, which commute with $H_{0}$. There are the relations [76]

$\displaystyle P_{0}\Omega=P_{0},\ \ \ \Omega P_{0}=\Omega.$ (S3)

Now we consider the perturbative expansion. The square root is expanded as
[95]

$\displaystyle(\Omega^{\dagger}\Omega)^{-1/2}$
$\displaystyle=P_{0}+\sum_{j=1}^{\infty}\frac{(-1)^{j}}{2^{2j}}\,{}_{2j}\textrm{C}_{j}\,(\Omega^{\dagger}\Omega-
P_{0})^{j}$ (S4) $\displaystyle\simeq
P_{0}-\frac{1}{2}\Omega_{1}^{\dagger}\Omega_{1},$ (S5)

where $\Omega$ is expanded as $\Omega=\Omega_{0}+\Omega_{1}+\cdots$ and only
the contributions up to second-order are kept. The second-order effective
Hamiltonian is given by

$\displaystyle\mathscr{H}_{\rm eff}$
$\displaystyle=P_{0}(H_{0}+V)P_{0}+\frac{1}{2}(P_{0}V\Omega_{1}+\Omega_{1}^{\dagger}VP_{0}).$
(S6)

The matrix element is evaluated as

$\displaystyle\langle a|\mathscr{H}_{\rm eff}|b\rangle=\langle
a|(H_{0}+V)|b\rangle$ $\displaystyle\hskip 14.22636pt+\frac{1}{2}\langle
a|V\Big{(}Q_{0}\frac{1}{E_{a}-H_{0}}+\frac{1}{E_{b}-H_{0}}Q_{0}\Big{)}V|b\rangle,$
(S7)

where $|a,b\rangle$ belong to the model Hilbert space.

## SM 3: SU($N$) generators

When we expand the effective Hamiltonian, we take the matrix representation
$O_{\alpha\beta}^{\xi}$ as SU($N$) generators, where
$\alpha,\beta,\xi\in\mathbb{N}$, $\alpha,\beta\in[1,N]$, and
$\xi\in[0,N^{2}-1]$. Just for convenience, we impose the Hermiticity,
completeness, and orthonormality for the matrix basis:

$\displaystyle(\hat{O}^{\xi})^{\dagger}=\hat{O}^{\xi},$ (S8)
$\displaystyle\sum_{\xi}(O^{\xi}_{\alpha\beta})^{*}O^{\xi}_{\alpha^{\prime}\beta^{\prime}}=\delta_{\alpha\alpha^{\prime}}\delta_{\beta\beta^{\prime}},$
(S9)
$\displaystyle\sum_{\alpha\beta}(O^{\xi}_{\alpha\beta})^{*}O^{\xi^{\prime}}_{\alpha\beta}=\Tr\hat{O}^{\xi}\hat{O}^{\xi^{\prime}}=\delta_{\xi\xi^{\prime}},$
(S10)

where the hat ($\hat{\ }$) symbol represents a matrix with respect to the
index $\alpha$. Then the coupling constant satisfies
$I_{ij}^{\xi\xi^{\prime}}=I_{ji}^{\xi^{\prime}\xi}\in\mathbb{R}$. The matrix
representation of SU($N$) generators consists of $N$ diagonal matrices and
$N^{2}-N$ off-diagonal ones [96]. The diagonal components are explicitly
written as

$\displaystyle
O_{\alpha\beta}^{\xi=0}=\frac{1}{\sqrt{N}}\delta_{\alpha\beta},$ (S11)
$\displaystyle
O_{\alpha\beta}^{\eta}=\frac{1}{2\sqrt{\eta(\eta+1)}}\quantity(\sum_{\zeta=1}^{\eta}\delta_{\alpha\zeta}\delta_{\beta\zeta}-\eta\delta_{\alpha,\eta+1}\delta_{\beta,\eta+1}),$
(S12)

where $\eta\in[1,N-1]$. $\hat{O}^{0}$ is proportional to the identity matrix.
The concrete forms of off-diagonal matrices are constructed by putting
$1/\sqrt{2}$ or $-\mathrm{i}/\sqrt{2}$ at one element of the upper triangular
block. The lower triangular block are determined from Hermiticity.

## SM 4: Definition of the physical quantities

In this section, we express the spin, magnetic orbital and quadrupole
operators, whose expectation values are calculated in the main text. The spin
operator is defined by

$\displaystyle\bm{S}_{i}=\frac{1}{2}\sum_{\gamma\sigma\sigma^{\prime}}c_{i\sigma\gamma}^{\dagger}\bm{\sigma}_{\sigma\sigma^{\prime}}c_{i\sigma^{\prime}\gamma},$
(S13)

where $\bm{\sigma}$ is the Pauli matrix. As for the orbital dependent
quantities, we take the orbital basis $(|xy\rangle,|yz\rangle,|zx\rangle)$
with the local coordinate of the $t_{2g}$ electrons where $z$-direction is
along the local three-fold rotationally symmetric axis at each site. The
magnetic orbital is written as

$\displaystyle\bm{L}_{i}=\sum_{\gamma\gamma^{\prime}\sigma}c_{i\sigma\gamma}^{\dagger}\bm{\ell}_{\gamma\gamma^{\prime}}c_{i\sigma\gamma^{\prime}},$
(S14)

where

$\displaystyle\hat{\ell}^{x}=\begin{pmatrix}0&0&-\mathrm{i}\\\ 0&0&0\\\
\mathrm{i}&0&0\end{pmatrix},\hat{\ell}^{y}=\begin{pmatrix}0&\mathrm{i}&0\\\
-\mathrm{i}&0&0\\\ 0&0&0\end{pmatrix},\hat{\ell}^{z}=\begin{pmatrix}0&0&0\\\
0&0&\mathrm{i}\\\ 0&-\mathrm{i}&0\end{pmatrix}.$ (S15)

We have omitted the site index $i$ to make the notation simple. We also write
the quadrupole moment as

$\displaystyle
Q_{i}^{\eta}=\sum_{\gamma\gamma^{\prime}\sigma}c_{i\sigma\gamma}^{\dagger}q_{\gamma\gamma^{\prime}}^{\eta}c_{i\sigma\gamma^{\prime}},$
(S16)

where the matrix representation $\hat{q}^{\eta}$ is constructed from
$\hat{\ell}^{\mu}$ as

$\displaystyle\hat{q}^{x^{2}-y^{2}}=(\hat{\ell}^{x})^{2}-(\hat{\ell}^{y})^{2},$
(S17)
$\displaystyle\hat{q}^{z^{2}}=\frac{1}{\sqrt{3}}(2(\hat{\ell}^{z})^{2}-(\hat{\ell}^{x})^{2}-(\hat{\ell}^{y})^{2}),$
(S18)
$\displaystyle\hat{q}^{xy}=\hat{\ell}^{x}\hat{\ell}^{y}+\hat{\ell}^{y}\hat{\ell}^{x},$
(S19)
$\displaystyle\hat{q}^{yz}=\hat{\ell}^{y}\hat{\ell}^{z}+\hat{\ell}^{z}\hat{\ell}^{y},$
(S20)
$\displaystyle\hat{q}^{zx}=\hat{\ell}^{z}\hat{\ell}^{x}+\hat{\ell}^{x}\hat{\ell}^{z}.$
(S21)

Combining the above matrices, we define another multipoles. The electric
dipole moment is written as

$\displaystyle
G_{i}^{\mu}=\frac{1}{2}\sum_{\nu\lambda}\varepsilon_{\mu\nu\lambda}\sum_{\gamma\gamma^{\prime}\sigma\sigma^{\prime}}c_{i\gamma\sigma}^{\dagger}\ell_{\gamma\gamma^{\prime}}^{\nu}\sigma_{\sigma\sigma^{\prime}}^{\lambda}c_{i\gamma^{\prime}\sigma^{\prime}},$
(S22)

where $\varepsilon_{\mu\nu\lambda}$ is the completely anti-symmetric tensor.
At last, the magnetic octupole moment is written as

$\displaystyle
T_{i}^{\eta\mu}=\frac{1}{2}\sum_{\gamma\gamma^{\prime}\sigma\sigma^{\prime}}c_{i\gamma\sigma}^{\dagger}q_{\gamma\gamma^{\prime}}^{\eta}\sigma_{\sigma\sigma^{\prime}}^{\mu}c_{i\gamma\sigma}.$
(S23)

The numerical results of the orbital dependent quantities $L,Q,G$ and $T$
shown in the main text are rotated to the globally defined axes (see the inset
of Fig. 1 of the main text).

## SM 5: Details of the classical model

###  .1 Explicit form of the coherent state

We consider the coherent state [52]

$\displaystyle|\Omega\rangle$
$\displaystyle=\sum_{\alpha=1}^{N}c_{\alpha}(\Omega)|\alpha\rangle,$ (S24)
$\displaystyle c_{\alpha}(\Omega)$
$\displaystyle=\mathrm{e}^{\mathrm{i}\varphi_{\alpha-1}}\cos\xi_{\alpha}\prod_{\beta=1}^{\alpha-1}\sin\xi_{\beta},$
(S25)

where $\varphi_{0,\cdots,N-1}$ and $\xi_{1,\cdots,N}$ with $\varphi_{0}=0$ and
$\xi_{N}=0$. Thus we have $2(N-1)$ parameters. Here, we have omitted the site
index $i$. It follows that

$\displaystyle\langle\Omega|\Omega\rangle=1,$ (S26)
$\displaystyle\frac{N!}{\pi^{N-1}}\int\mathrm{d}\Omega\,|\Omega\rangle\langle\Omega|=1,$
(S27)
$\displaystyle\mathrm{d}\Omega=\prod_{\alpha=1}^{N-1}\cos\xi_{\alpha}\sin^{2(N-\alpha)-1}\xi_{\alpha}\
\mathrm{d}\xi_{\alpha}\mathrm{d}\varphi_{\alpha},$ (S28)
$\displaystyle\innerproduct{\Omega}{\partial\Omega}=\mathrm{i}\sum_{\alpha=1}^{N}\absolutevalue{c_{\alpha}(\Omega)}^{2}\partial\varphi_{\alpha}.$
(S29)

The parameters are $\xi_{1},\cdots,\xi_{N-1}\in[0,\pi/2]$ and
$\varphi_{1},\cdots,\varphi_{N-1}\in[0,2\pi)$ [52].

###  .2 Over-relaxation update

We provide a detailed explanation for the over-relaxation update. We begin
with the local Hamiltonian

$\displaystyle\mathcal{H}_{{\rm loc}}$
$\displaystyle=-\sum_{\xi}\tilde{H}^{\xi}\mathcal{O}(\Omega)$ (S30)
$\displaystyle=\sum_{\alpha\beta}h_{\alpha\beta}c_{\alpha}^{*}(\Omega)c_{\beta}(\Omega)=\sum_{\gamma}\Lambda_{\gamma}|d_{\gamma}(\Omega)|^{2},$
(S31)

where the site index is omitted for simplicity. The diagonalization is
performed in the right-most side by the unitary matrix $V$:

$\displaystyle
d_{\gamma}(\Omega)=\sum_{\alpha}V^{\dagger}_{\gamma\alpha}c_{\alpha}(\Omega).$
(S32)

It is apparent at this point that the energy does not change by the phase
transformation

$\displaystyle d_{\gamma}(\Omega)\to
d_{\gamma}\mathrm{e}^{\mathrm{i}\theta_{\gamma}}=\sum_{\alpha}\mathrm{e}^{\mathrm{i}\theta_{\gamma}}V^{\dagger}_{\gamma\alpha}c_{\alpha}(\Omega)\equiv\sum_{\alpha}V^{\dagger}_{\gamma\alpha}c_{\alpha}(\Omega^{\prime}),$
(S33)

with which the coherent state is transformed as $\Omega\to\Omega^{\prime}$.
For the efficient update to $\Omega^{\prime}$, we would like to know the
“opposite side” of $\Omega$ [85]. For this purpose, we minimize the inner
product defined by

$\displaystyle\langle\Omega|\Omega^{\prime}\rangle$
$\displaystyle=\sum_{\gamma}\mathrm{e}^{\mathrm{i}\theta_{\gamma}}|d_{\gamma}(\Omega)|^{2}.$
(S34)

We define the norm

$\displaystyle\mathcal{N}[\theta]$
$\displaystyle=|\langle\Omega|\Omega^{\prime}\rangle|^{2},$ (S35)

and we find that
$\displaystyle\frac{\partial\mathcal{N}}{\partial\theta_{\gamma}}=0$ is
satisfied if $\mathrm{e}^{\mathrm{i}\theta_{\gamma}}=\pm 1\equiv s_{\gamma}$.
The set of signs ($s_{1},\cdots,s_{N}$) is determined so as to minimize the
norm $\mathcal{N}$. We search for the solution by considering $2^{N}$
possibilities, which is same as the bipartitioning problem. This procedure
reproduces the over-relaxation update usually used for the SU(2) case.

##  SM 6: Results for simple SU(2) model

Since the model with general $N$ is complicated, it is useful to summarize the
results for the $N=2$ case as a benchmark. Here we consider the single-orbital
Hubbard model on the cubic lattice. The Hamiltonian is given by

$\displaystyle\mathscr{H}$ $\displaystyle=-t\sum_{\langle
ij\rangle\sigma}c_{i\sigma}^{\dagger}c_{j\sigma}+U\sum_{i}n_{i\uparrow}n_{i\downarrow}.$
(S36)

Applying the second-order perturbation theory in the strong coupling limit, we
obtain the effective Hamiltonian for $n=1$ as

$\displaystyle\mathscr{H}_{\rm eff}$ $\displaystyle=I\sum_{\langle
ij\rangle}\sum_{\xi=1}^{3}\mathscr{O}^{\xi}_{i}\mathscr{O}^{\xi}_{j},$ (S37)

where $I=2t^{2}/U$ and

$\displaystyle\mathscr{O}^{\xi}=\sum_{\alpha,\beta=1}^{2}|\alpha\rangle_{i}\,O^{\xi}_{\alpha\beta}\,{}_{i}\langle\beta|.$
(S38)

The state vectors in model space are
$|1\rangle_{i}=c_{i\uparrow}^{\dagger}|0\rangle_{i}$ and
$|2\rangle_{i}=c_{i\downarrow}^{\dagger}|0\rangle_{i}$, where $|0\rangle_{i}$
is the vacuum at site $i$. The matrices are constructed following the
procedure in SM3:

$\displaystyle O^{1}=\frac{1}{\sqrt{2}}\begin{pmatrix}1&0\\\
0&-1\end{pmatrix},O^{2}=\frac{1}{\sqrt{2}}\begin{pmatrix}0&1\\\
1&0\end{pmatrix},O^{3}=\frac{1}{\sqrt{2}}\begin{pmatrix}0&-\mathrm{i}\\\
\mathrm{i}&0\end{pmatrix}.$ (S39)

We note that the above procedure is performed automatically in our framework.

Figure S1:  (a) Electronic band structure and (b) spin excitation spectrum
${\rm Im\,}\chi(\bm{q},\omega)/\omega$ for the Heisenberg model on a cubic
lattice. The parameters are chosen as $t=1$ eV, $U=8$ eV and $T=10^{-3}$ eV.

We take the two nearest neighbor atoms in the unit cell, which are labeled by
the sublattice index A,B. In this case, the antiferromagnetism occurs as a
$\bm{q}=\bm{0}$ state. Here, we solve the Heisenberg model by the mean-field
theory. The spin excitation spectra are shown in Fig. S1, which is an analog
of Fig. 2 of the main text. The gapless magnon mode around $\Gamma$
($\bm{q}=\bm{0}$) is clearly seen as expected.

The coherent state coefficients can also be explicitly written down
($\Omega_{1}=\xi\in[0,\pi/2]$, $\Omega_{2}=\varphi\in[0,2\pi)$):

$\displaystyle c_{1}(\Omega)$ $\displaystyle=\cos\xi,$ (S40) $\displaystyle
c_{2}(\Omega)$ $\displaystyle=\mathrm{e}^{\mathrm{i}\varphi}\sin\xi.$ (S41)

For the SU(2) case, $\varphi$ is interpreted as an azimuthal angle on the
Bloch sphere, and $\theta=2\xi\in[0,\pi]$ as a polar angle. The Berry
curvature matrix is given by

$\displaystyle\mathcal{B}$ $\displaystyle=\begin{pmatrix}0&\mathrm{i}\sin
2\xi\\\ -\mathrm{i}\sin 2\xi&0\end{pmatrix},$ (S42)

which leads to the Bloch’s equation of motion for spin dynamics.

References
See the list of references in the main text.