# Electronic origin of ferroic quadrupole moment under antiferroic quadrupole
orders and finite magnetic moment in $J_{\rm eff}=3/2$ systems

Haruhiro Kubo  Takayuki Ishitobi  Kazumasa Hattori  Department of Physics,
Tokyo Metropolitan University,
1-1, Minami-osawa, Hachioji, Tokyo 192-0397, Japan

###### Abstract

We study the electronic origin of parasitic ferroic quadrupole moments in
antiferroic quadrupole orders by extending a model studied in G. Chen et al.,
Phys. Rev. B 82, 174440 (2010) with the effective angular momentum $J_{\rm
eff}=3/2$ quartet ground states. Taking into account the first crystalline-
electric-field (CEF) excited doublet, cubic anisotropy in the quadrupole
moments emerges, which leads to the induced ferroic quadrupole moments in the
antiferro quadrupolar phases. The hybridization with the CEF excited quartet
states also causes finite magnetic moments compatible to the observed size of
the effective moment in typical $J_{\rm eff}=3/2$ systems, as opposed to the
naive expectation of vanishing moments in the $J_{\rm eff}=3/2$ systems. These
results suggest the importance of the corrections arising from the high-energy
CEF excited states in the $J_{\rm eff}=3/2$ systems.

###### pacs:

Valid PACS appear here

##  I Introduction

Strongly correlated electrons systems with orbital degrees of freedom have
attracted great attention in recent years Tokura and Nagaosa (2000). The
orbital degrees of freedom possess potential functions alternative to the
modern magnetic-based devicesFiebig  _et al._ (2016). Combining the orbital
and conventional spin degrees of freedom inevitably leads to the notion of
multipole momentsKuramoto  _et al._ (2009). Such multipole moments play
important role in correlated electron systems with strong spin-orbit
couplings. They can exhibit various fascinating phenomena such as topological
spin-orbital Mott insulatorsPesin and Balents (2010), spin-orbital liquid
statesJackeli and Khaliullin (2009); Nasu  _et al._ (2014); Corboz  _et al._
(2012); Kitagawa  _et al._ (2018), superconductivity mediated by multipolar
spin-orbital fluctuationsNomoto  _et al._ (2016); Hattori  _et al._ (2017),
and in doped Mott insulatorsWatanabe  _et al._ (2013).

For activating multipole physics, one needs high symmetry and strong spin-
orbit couplings. Possible candidates are Mott insulators with $d^{1}$
configuration surrounded by a regular oxygen octahedron. Although materials
with $3d$ electrons often exhibit the Jahn-Teller distortion as lowering
temperature, some 5$d$ electron systems keep their cubic symmetry even at low
temperatures. As a such $d$ electron system, double-perovskite compounds
$A_{2}B$$Tr$O6 ($A$=Ba, Sr,Ca; $B$= Mg, Ca, Sr, Ba, Zn, Cd, $Tr=$Re, Os,
Mo)Vasala and Karppinen (2015); Longo and Ward (1961); Wiebe  _et al._ (2003);
Sleight  _et al._ (1962); Chamberland and Levasseur (1979); Abakumov  _et al._
(1997); Wiebe  _et al._ (2002, 2003); Bramnik  _et al._ (2003); Yamamura  _et
al._ (2006); Marjerrison  _et al._ (2016); Hirai and Hiroi (2019) and Ta
chlorides $A_{2}$TaCl6 ($A=$K, Rb, Cs)Ishikawa  _et al._ (2019) have been
recently studied intensively. The oxygen octahedron crystalline-electric-field
(CEF) lifts the ten-fold degeneracy of the $d$ electrons to the high-energy
$e_{g}$ and the low-energy $t_{2g}$ states. The latter is split further by the
spin-orbit coupling and form so-called effective angular momentum $J_{\rm
eff}=3/2$ ground state and $J_{\rm eff}=1/2$ first excited state. In these
compounds, the spin-orbit coupling $\lambda$ is of the order of $\sim$0.2–0.4
eV, while the CEF gap between the ground $t_{2g}$ and the excited $e_{g}$
states $D$ is $D\sim$3–5 eV Paramekanti  _et al._ (2018); Oh  _et al._ (2018);
Ishikawa  _et al._ (2019).

What is unique to these systems is that the ground state with $J_{\rm
eff}=3/2$ has no magnetic dipole moment as a result of exact cancellation
between the spin and orbital angular momenta. This leads to various
possibilities of multipole orders at low temperatures. For example, Re based
double-perovskites exhibit ferro and antiferro quadrupole orders in addition
to dipole-octupole magnetic orders. In Ba2$Tr$OsO6 ($Tr=$Zn, Mg, and Ca) with
$d^{2}$ configuration, octupole orders have been suggested
recently.Paramekanti  _et al._ (2020); Maharaj  _et al._ (2020)

The vanishing magnetic dipole moment is also reflected in the small moment
typically $\sim 0.3$–$0.8\mu_{\rm B}$ in their high-temperature magnetic
susceptibility, where $\mu_{\rm B}$ is the Bohr magneton. This small but
finite value has been interpreted as a partial cancellation of the spin and
orbital angular momenta owing to the hybridization between the $d$ and the
oxygen $p$ orbitals.Ahn  _et al._ (2017); Hirai and Hiroi (2019); Ishikawa
_et al._ (2019)

In the pioneering work by Chen et al.,Chen  _et al._ (2010) a model containing
the ground-state CEF quartet originating from the $t_{2g}$ orbital was
constructed in the limit of the strong spin-orbit coupling $\lambda\to\infty$.
The quartet can be described by the effective angular momentum $J_{\rm
eff}=3/2$ and they predicted several interesting multipolar ordered states in
terms of the multipole moments of the $J_{\rm eff}=3/2$ multiplet: an
antiferroic quadrupole order with the $e_{g}$ irreducible representation and
ferro- and antiferro-magnetic dipole-octupole orders. These results are indeed
supported by the experimental observation of the phase transitions e.g., in
Ba2MgReO6Hirai and Hiroi (2019) and Ba2NaOsO6.Lu  _et al._ (2017)

Recently, Hirai et al., reported that ferroic quadrupole moments $\sim
O_{20}=3z^{2}-r^{2}$ emerge under the antiferroic order of the type
$O_{22}=x^{2}-y^{2}$ below $T_{q}=33$ K by their x-ray experimentsHirai  _et
al._ (2020). The presence of the ferroic moment can be understood by a simple
symmetry argument; the free energy for the $e_{g}$ orbital moments contains a
cubic coupling $\sim O_{22}^{2}O_{20}$. This induces ferroic quadrupole
moments proportional to the square of the antiferroic ones. They discuss this
can be due to the Jahn-Teller effects and the lattice anharmonicity. The
former has been recently analyzed and successfully reproduced the emergence of
the ferroic quadrupole momentsIwahara and Chibotaru (2022).

In this paper, we point out that $d$ electron CEF excited states can influence
the CEF ground quartet with $J_{\rm eff}=3/2$, which has no magnetic dipole
moment and no cubic coupling $\sim O_{22}^{2}O_{20}$ without couplings to
other degrees of freedom. In the $J_{\rm eff}=3/2$ multiplets, the local cubic
anisotropy which arises from the CEF potential vanishes. The anisotropy can
emerge when the CEF excited states are taken into account. Such anisotropy due
to the excited spin-singlet state ($\Gamma_{1}$ in the cubic symmetry) are
indeed discussed in the $\Gamma_{3}$ non-Kramers doublet ground state in Pr-
based compoundsHattori and Tsunetsugu (2014, 2016); Ishitobi and Hattori
(2021); Tsunetsugu  _et al._ (2021); Hattori  _et al._ (2022) and the analysis
there is also applicable to the case of the $d^{1}$ systems. This is because
the $J_{\rm eff}=3/2$ states is classified as $\Gamma_{8}$ irreducible
representation (irrep) in the cubic symmetry and can be regarded as a product
of the spin-1/2 and the orbital $e_{g}(\Gamma_{3})$. Indeed, the CEF first-
excited state is the $\Gamma_{7}$ state which is the spin-1/2 and the orbital-
singlet state. Thus, apart from the spin degrees of freedom, the orbital
sector is identical in the two systems. For the small but finite magnetic
dipole moment, it will be shown that the corrections of order $\lambda/D\sim
0.1$ leads to non-negligible contribution to the magnetic moment in realistic
systems in this paper.

This paper is organized as follows. In Sec. II, we review the model proposed
in Ref. Chen  _et al._ , 2010 and explicitly introduce the matrix form of
quadrupole operators including the excited states. The exchange Hamiltonian is
rewritten in terms of these quadrupole and spin-orbital operators to make this
paper self-contained form. In Sec. III, we show the results of the two-site
mean-field approximation. We discuss the effects of the excited CEF state on
the ferro components of the order parameter in the antiferroic quadrupole
ordered state and also on the phase-diagram. The temperature-magnetic field
phase diagrams are also analyzed. In Sec. IV, we discuss the results in this
paper and related materials. We also discuss how the finite magnetic moment
emerges in the ground state $J_{\rm eff}=3/2$ state. Finally, Sec. V
summarizes this paper.

##  II Model

We start by introducing the general model describing $d^{1}$ electron
configuration for arbitrary spin-orbit coupling and CEF strength. This means
that both $t_{2g}$ and $e_{g}$ degrees of freedom together with the spin $1/2$
ones are taken into account. We then derive an effective $t_{2g}$-dominant
model with six states, ignoring the excited $e_{g}$ dominant states. The
interaction between the effective $t_{2g}$ electrons are introduced as
similarly to the study in Ref. Chen  _et al._ , 2010. We will not project them
onto the effective total angular momentum $J_{\rm eff}=3/2$, but keep both
$J_{\rm eff}=3/2$ and $1/2$ constructed by the $t_{2g}$ orbitals.

###  II.1 Local Hamiltonian

The local part of the Hamiltonian is a conventional one with the CEF parameter
$B_{4}$ and the spin-orbit coupling $\lambda>0$ as

$\displaystyle H_{\rm loc}$ $\displaystyle=H_{\rm CEF}+H_{\rm SO},$ (1)
$\displaystyle H_{\rm CEF}$
$\displaystyle=\frac{B_{4}}{2}\sum_{i,\sigma,m,m^{\prime}}\gamma_{mm^{\prime}}d_{m\sigma}^{\dagger}d_{m^{\prime}\sigma},$
(2) $\displaystyle H_{\rm SO}$
$\displaystyle=\lambda\sum_{i,\sigma,m,m^{\prime}}({\bm{L}})_{mm^{\prime}}\cdot\bm{s}_{\sigma\sigma^{\prime}}d_{m\sigma}^{\dagger}d_{m^{\prime}\sigma^{\prime}},$
(3)

Here, $d_{m\sigma}$ is the $d$ electron annihilation operator with the $z$
component of the orbital angular momentum $m=0,\pm 1,\pm 2$ and the spin
$\sigma=\uparrow,\downarrow$. The CEF potential for the $d$ electrons are
parameterized by
$\gamma_{mm^{\prime}}\equiv(\frac{15}{2}m^{2}-\frac{35}{2}|m|+6)\delta_{m,m^{\prime}}+5\delta_{|m-m^{\prime}|,4}$,
where $\delta_{m,m^{\prime}}$ is the Kronecker delta. $\bm{L}$ and $\bm{s}$
are the orbital angular momentum and the spin-1/2 matrices for the $d$
electrons, respectively.

The eigenstates of $H_{\rm loc}$ are split into $4+2+4$ and their eigenvalues
are given as

$\displaystyle\epsilon_{1/2}$ $\displaystyle=-2B_{4}+\lambda,$ (4)
$\displaystyle\epsilon_{3/2}$
$\displaystyle=\frac{1}{4}\left(2B_{4}-\lambda-\sqrt{5}\sqrt{20B_{4}^{2}+4B_{4}\lambda+5\lambda^{2}}\right),$
(5) $\displaystyle\epsilon_{3/2}^{\prime}$
$\displaystyle=\frac{1}{4}\left(2B_{4}-\lambda+\sqrt{5}\sqrt{20B_{4}^{2}+4B_{4}\lambda+5\lambda^{2}}\right),$
(6)

where the subscripts represent the effective angular momenta; $3/2$ ($1/2$)
corresponds to $\Gamma_{8}$($\Gamma_{7}$) state in the cubic point group Oh.

For small $\lambda/B_{4}$ with $B_{4}>0$,

$\displaystyle\epsilon_{3/2}$
$\displaystyle\simeq-2B_{4}-\frac{\lambda}{2}-\frac{3\lambda^{2}}{10B_{4}}+\cdots,$
(7) $\displaystyle\epsilon_{3/2}^{\prime}$ $\displaystyle\simeq
3B_{4}+\frac{3\lambda^{2}}{10B_{4}}+\cdots.$ (8)

When the CEF potential is sufficiently large, $\epsilon_{3/2}^{\prime}$ is
much larger than the other two and the excitation gap $D$ is $D\simeq
5B_{4}\sim 3$–$5$ eV for the compounds mentioned in the
Introduction.Paramekanti  _et al._ (2018); Oh  _et al._ (2018); Ishikawa  _et
al._ (2019) Thus, the quartet with its eigenenergy $\epsilon_{3/2}$
corresponds to the ground state with $J_{\rm eff}=3/2$ and the first excited
state is that with $\epsilon_{1/2}$. The former approximately consists of the
$t_{2g}$ electrons, while the latter is purely $t_{2g}$ origin, see Appendix
A. In the following, we will concentrate on these six states and ignore the
higher energy states at $\epsilon_{3/2}^{\prime}$. The matrices which will be
shown in the following [Eqs. (9),(11), and (12)] are calculated within the
full ten-dimensional Hilbert space with keeping $O(\varepsilon^{2})$ terms for
later purposes. The six states can be labeled by the diagonal elements of the
total angular momentum $J_{z}$ for small $\varepsilon\equiv\lambda/B_{4}\sim
0.5$:

$\displaystyle J_{z}$
$\displaystyle\simeq\frac{1}{3}\begin{bmatrix}\frac{5}{2}&0&-2\sqrt{2}+\alpha_{J}&0&0&0\\\
&-\frac{5}{2}&0&-2\sqrt{2}+\alpha_{J}&0&0\\\ &&\frac{1}{2}+\beta_{J}&0&0&0\\\
&&0&-\frac{1}{2}-\beta_{J}&0&0\\\ &&&&-\frac{3}{2}&0\\\ &&&&&\frac{3}{2}\\\
\end{bmatrix}.$ (9)

Here the empty parts in the matrix have been omitted since
$J_{z}=J_{z}^{\dagger}$. The constants $\alpha_{J}$ and $\beta_{J}$ are given
as
$\alpha_{J}=\frac{3\sqrt{2}}{5}\varepsilon+O\left(\varepsilon^{3}\right),\quad\beta_{J}=\frac{12}{5}\varepsilon-\frac{9}{25}\varepsilon^{2}+O\left(\varepsilon^{3}\right).$
Note that the factor $1/3$ in Eq. (9) and the basis of the matrix $J_{z}$ is
the eigenstates for $H_{\rm loc}={\rm
diag}(\epsilon_{1/2},\epsilon_{1/2},\epsilon_{3/2},\epsilon_{3/2},\epsilon_{3/2},\epsilon_{3/2})$.
See the wavefunctions in Appendix A. From this expression, it is clear that
the quartet can be described by the effective angular momentum $J_{\rm
eff}=3/2$ with $J_{\rm eff}^{z}=3J_{z}$ when the excited states are ignored
for $\lambda\to\infty$ with $\varepsilon=\lambda/B_{4}\to 0$.

For later purposes, it is useful to show the operators constructed by the
$t_{2g}$ occupation number $n_{xy},\ n_{yz}$, and $n_{zx}$; $n\equiv
n_{xy}+n_{yz}+n_{zx}$, $u\equiv 2n_{xy}-n_{yz}-n_{zx}$, and
$v\equiv\sqrt{3}(n_{yz}-n_{zx})$.

$\displaystyle n$ $\displaystyle\simeq{\rm
diag}(1,1,1-\delta,1-\delta,1-\delta,1-\delta),$ (10) $\displaystyle u$
$\displaystyle\simeq\small\begin{bmatrix}0&0&\frac{2-\delta}{\sqrt{2}}&0&0&0\\\
&0&0&-\frac{2-\delta}{\sqrt{2}}&0&0\\\ &&1-\delta&0&0&0\\\ &&&1-\delta&0&0\\\
&&&&-1+\delta&0\\\ &&&&&-1+\delta\\\ \end{bmatrix},$ (11) $\displaystyle v$
$\displaystyle\simeq\small\begin{bmatrix}0&0&0&0&\frac{2-\delta}{\sqrt{2}}&0\\\
&0&0&0&0&\frac{2-\delta}{\sqrt{2}}\\\ &&0&0&-1+\delta&0\\\ &&&0&0&1-\delta\\\
&&&&0&0\\\ &&&&&0\\\ \end{bmatrix},$ (12)

where $\delta=\frac{3}{50}\varepsilon^{2}+O(\varepsilon^{3})$, and as
similarly to Eq. (9), we have omitted the lower triangle part. For $\delta\to
0$, $n\to\mathbbm{1}$ (an identity matrix) and $e_{g}$ quadrupole moments
$(u,v)$ are reduced to the quadrupole moments $(q_{u},q_{v})$ constructed by
the angular momentum ${\bm{L}}$ as

$\displaystyle q_{u}=\frac{2L_{z}^{2}-L_{x}^{2}-L_{y}^{2}}{3},\quad
q_{v}=\frac{\sqrt{3}(L_{x}^{2}-L_{y}^{2})}{3}.$ (13)

Note that the $u$ and $v$ include the virtual processes accross the CEF
excited quartet states, which is different from that calculated by the
$L_{x,y,z}$ restricted within the $t_{2g}$ states. It is natural to obtain
these expressions in terms of $\bm{L}$ for $\delta\to 0$ instead of those of
the total angular momentum $\bm{J}$ since $n_{xy,yz,zx}$ does not depend on
the spin. We note that $(u,v)$ possess the offdiagonal matrix elements
$\pm\frac{2-\delta}{\sqrt{2}}$ between the ground states and the excited
states and their magnitudes are larger than those among the ground states.
This is the source of the anisotropy in the quadrupole moments and enables one
to obtain ferro quadrupole moments under antiferro quadrupole orders. In Fig.
1, the full expressions, i.e., without assumption of the small
$\varepsilon(=\lambda/B_{4})$ are shown as a function of $\varepsilon$. The
offdiagonal elements for $u$ and $q_{u}$ are exactly the same, while the
diagonal ones differ as $\varepsilon$ increases. For the realistic parameter
regime $\varepsilon\lesssim 0.5$, the difference between the full expression
and the approximated one shown in Eqs. (11) and (12) are quantitatively the
same.

The finite offdiagonal elements $\pm(2-\delta)/\sqrt{2}\sim\pm\sqrt{2}$ for
$\varepsilon\lesssim 0.5$ leads to a finite cubic anisotropy in the local
quadrupole free energy, which can be calculated by the local CEF model via the
Legendre transformationIshitobi and Hattori (2023) as

$\displaystyle F_{q}^{\rm loc}\sim F_{q0}^{\rm
loc}+a(\phi_{u}^{2}+\phi_{v}^{2})-b(\phi_{u}^{3}-3\phi_{u}\phi_{v}^{2})+\cdots,$
(14)

where $F_{q0}^{\rm loc},a$, and $b$ are constants and $\phi_{u}$ and
$\phi_{v}$ correspond to the quadrupole fields for $u$ and $v$, respectively.
The coefficient of the cubic anisotropy $b$ is given by

$\displaystyle
b\simeq\frac{1}{{2}\beta^{2}(\frac{3}{2}\lambda)}\left(\frac{2-\delta}{\sqrt{2}}\right)^{2}.$
(15)

Here, $\beta$ is the inverse of temperature $T$. This is the consequence of
standard Landau expansion of quadrupole free energy for $\beta\lambda\gg 1$.
The $T^{2}$ dependence of $b$ is meaningful near the quadrupolar transition
temperature $T_{q}\simeq 33$ K and $T$ is usually replaced by $T_{q}$ for its
phenomenological analysis. We note that $b$ is proportional to the square of
the offdiagonal element in $u$ and $v$ in Eqs. (11) and (12). The denominator
$\frac{3}{2}\lambda$ represents that the perturbative processes to the CEF
excited doublet are important. As for the opposite limit $\beta\lambda\sim 0$,
$b\propto T$.

Figure 1: The matrix elements of $u$ and $q_{u}$ as a function of
$\varepsilon=\lambda/B_{4}$. For the offdiagonal elements, $u$ and $q_{u}$ are
identical: $(u)_{13}=-(u)_{24}=(q_{u})_{13}=-(q_{u})_{24}$, while for the
diagonal ones, they are different for finite $\varepsilon$:
$(u)_{33}=(u)_{44}=-(u)_{55}=-(u)_{66}$.
$(q_{u})_{33}=(q_{u})_{44}=-(q_{u})_{55}=-(q_{u})_{66}$. The approximated
values for the matrix elements of $u$: $1-\delta$ and $(2-\delta)/\sqrt{2}$,
are also plotted.

###  II.2 Interactions in $t_{2g}$ manifolds

In this subsection, we introduce the exchange interactions between the
$t_{2g}$ electrons. The model used in this paper is basically given in Ref.
Chen  _et al._ , 2010, but for making this paper self-contained we summarize
the model in the following.

First, we discuss exchange interactions in the quadrupole sector. We use the
quadrupole-quadrupole interactions introduced in Ref. Chen  _et al._ , 2010
with slightly different notationTsunetsugu  _et al._ (2021); Hattori  _et al._
(2022) including two parameters $g_{\rm iso}$ and $g_{\rm ani}$ as

$\displaystyle H_{q}=\sum_{\langle i,j\rangle}(g_{\rm
iso}\bm{Q}_{i}\cdot\bm{Q}_{j}+g_{\rm
ani}\bm{Q}_{i}\cdot{\mathsf{K}}_{ij}\bm{Q}_{j}).$ (16)

Here, ${\bm{Q}}_{i}\equiv(u_{i},v_{i})$ is the $e_{g}$ quadrupole vector at
the $i$th site and the sum runs over the nearest-neighbor sites on the fcc
lattice. The bond directional anisotropic interactions are parameterized by
the matrix ${\mathsf{K}}_{ij}$:

$\displaystyle{\mathsf{K}}_{ij}$ $\displaystyle=\begin{bmatrix}-1&0\\\
0&1\end{bmatrix}$ $\displaystyle\equiv{\mathsf{K}}^{3}$ $\displaystyle\quad
ij\ \text{bond}\parallel\text{$xy$ plane},$ (17)
$\displaystyle{\mathsf{K}}_{ij}$ $\displaystyle=\begin{bmatrix}-c&s\\\
s&c\end{bmatrix}$ $\displaystyle\equiv{\mathsf{K}}^{1}$ $\displaystyle\quad
ij\ \text{bond}\parallel\text{$yz$ plane},$ (18)
$\displaystyle{\mathsf{K}}_{ij}$ $\displaystyle=\begin{bmatrix}-c&-s\\\
-s&c\end{bmatrix}$ $\displaystyle\equiv{\mathsf{K}}^{2}$ $\displaystyle\quad
ij\ \text{bond}\parallel\text{$zx$ plane},$ (19)

with $c\equiv\cos\frac{2\pi}{3}$ and $s\equiv\sin\frac{2\pi}{3}$. One can also
rewrite them as

${\mathsf{K}}^{n}=-\hat{t}_{n}\hat{t}_{n}^{\ \rm T},\quad\hat{t}_{n}^{\ \rm
T}\equiv(\cos n\omega,\sin n\omega),\quad n=1,2,3,$ (20)

with $\omega\equiv 2\pi/3$. Note that the unit vector $\hat{t}_{1,2,3}$
represents the projection operator to $3x^{2}-r^{2}$, $3y^{2}-r^{2}$, and
$3z^{2}-r^{2}$ type orbital, respectively. The parameter $V$ in Ref. Chen  _et
al._ , 2010 corresponds to $g_{\rm iso}=-7V/72$ and $g_{\rm ani}=-25V/72$.
Since the ratio $g_{\rm ani}/g_{\rm iso}=25/7>2$, this naively suggests an
antiferro “$v$” order (O22 type) with the ordering vector at the X point:
$(0,0,2\pi)$ from the results for $E_{g}$ non-Kramers doublet systems on the
fcc lattice.Tsunetsugu  _et al._ (2021); Hattori  _et al._ (2022) For a more
general situation, the ratio $g_{\rm ani}/g_{\rm iso}\neq 25/7$, but in this
paper we restrict ourselves to the case derived in Ref. Chen  _et al._ , 2010
since the modification in $g_{\rm ani}/g_{\rm iso}$ leads to just a
qualitative difference for our purpose in this paper.

For the spin part, by using the spin operator for the $\rho=xy,yz$, and $zx$
orbital:
$S_{i,\rho}^{\mu}\equiv\frac{1}{2}\sum_{\sigma\sigma^{\prime}}d^{\dagger}_{i,\rho\sigma}(\hat{\sigma}^{\mu})_{\sigma\sigma^{\prime}}d_{i,\rho\sigma^{\prime}}$,
antiferromagnetic interactions are given as

$\displaystyle H_{s1}=$ $\displaystyle J\sum_{\langle
i,j\rangle}^{xy\text{-plane}}\left(\sum_{\mu}S_{i,xy}^{\mu}S_{j,xy}^{\mu}-\frac{1}{4}n_{i,xy}n_{j,xy}\right)$
$\displaystyle+\left(xy\ \to\ yz\ {\rm and}\ zx\right).$ (21)

Here,
$n_{i,\rho}\equiv\sum_{\sigma}d^{\dagger}_{i,\rho\sigma}d_{i,\rho\sigma}$ is
the number operator for $\rho=xy,yz$, and $zx$ orbitals. There are also
ferromagnetic interaction,

$\displaystyle H_{s2}=$ $\displaystyle-J^{\prime}\\!\\!\\!\sum_{\langle
i,j\rangle}^{xy\text{-plane}}\left[\sum_{\mu}S_{i,xy}^{\mu}(S_{j,yz}^{\mu}+S_{j,zx}^{\mu})+i\leftrightarrow
j\right]$ $\displaystyle+\frac{3J^{\prime}}{2}\sum_{\langle i,j\rangle\in
xy\text{-plane}}n_{i,xy}n_{j,xy}$ $\displaystyle+\left(xy\ \to\ yz\ {\rm and}\
zx\right).$ (22)

$n_{i\rho}$ is represented by the quadrupole operators as

$\displaystyle n_{i,xy}$ $\displaystyle=\frac{1}{3}(n_{i}+u_{i}),$ (23)
$\displaystyle n_{i,yz}$
$\displaystyle=\frac{1}{3}\left(n_{i}-\frac{1}{2}u_{i}+\frac{\sqrt{3}}{2}v_{i}\right),$
(24) $\displaystyle n_{i,zx}$
$\displaystyle=\frac{1}{3}\left(n_{i}-\frac{1}{2}u_{i}-\frac{\sqrt{3}}{2}v_{i}\right),$
(25)

The terms consisting of $n_{i,\rho}$’s in Eqs. (21) and (22) can be rewritten
by the quadrupole forms as

$\displaystyle H^{\prime}_{q}=-\frac{J-6J^{\prime}}{72}\sum_{\langle
i,j\rangle}\bm{Q}_{i}\cdot(\mathbbm{1}-\mathsf{K}_{ij})\bm{Q}_{j},$ (26)

which renormalizes $g_{\rm iso}$ and $g_{\rm ani}$ in Eq. (16).

Summing up Eqs. (16), (21), and (22) with Eq. (26), we obtain the total
nearest-neighbor exchange Hamiltonian as

$\displaystyle H_{\rm int}=$ $\displaystyle\sum_{\langle
i,j\rangle}\left(\tilde{g}_{\rm iso}\bm{Q}_{i}\cdot\bm{Q}_{j}+\tilde{g}_{\rm
ani}\bm{Q}_{i}\cdot{\mathsf{K}}_{ij}\bm{Q}_{j}\right)$
$\displaystyle+(J+2J^{\prime})\sum_{\rho}\sum_{\langle
i,j\rangle}^{\rho{\text{-plane}}}\bm{S}_{i,\rho}\cdot\bm{S}_{j,\rho}$
$\displaystyle-J^{\prime}\sum_{\rho}\sum_{\langle
i,j\rangle}^{\rho\text{-plane}}(\bm{S}_{i,\rho}\cdot\bm{S}_{j}+\bm{S}_{j,\rho}\cdot\bm{S}_{i}),$
(27)

where $\rho=xy,yz$, and $zx$. We have introduced the $t_{2g}$ spin operators
at the $i$ site: $\bm{S}_{i}=\sum_{\rho}\bm{S}_{i,\rho}$ and the renormalized
quadrupole interactions $\tilde{g}_{\rm iso}=g_{\rm
iso}-\frac{1}{72}(J-6J^{\prime})$ and $\tilde{g}_{\rm ani}=g_{\rm
ani}+\frac{1}{72}(J-6J^{\prime})$.

##  III Mean-field results

In this section, we show the results of the mean-field analysis of the model
(27) with the local part of the Hamiltonian (1). In Sec. III.1, we demonstrate
the results for the magnetic field $\bm{h}={\bf 0}$. For ${\bm{h}}={\bf 0}$,
the CEF excited quartet plays just a minor role as is evident from the small
factors $\delta\simeq\frac{3}{50}\varepsilon^{2}$ and
$\alpha\simeq\frac{3\sqrt{2}}{100}\varepsilon^{2}$ in the matrix elements in
the Hamiltonian (27) in addition to the energy gap $D\sim 5B_{4}$. See Eqs.
(11), (12), and (56)–(58). Thus, one can consider the
$\varepsilon=\lambda/B_{4}\to 0$ limit neglecting the CEF excited quartet
states with setting $\epsilon_{3/2}=0$, $\epsilon_{1/2}=\frac{3\lambda}{2}$,
and $\epsilon_{3/2}^{\prime}\simeq 5B_{4}\to\infty$ in Eqs. (4), (7), and (8),
respectively. Here, we shift the energy in order to set the energy of the CEF
ground quartet states to 0. The matrix elements of various operators are also
simplified in this limit: These simplifications are justified for
${\bm{h}}={\bf 0}$. In Sec. III.2, we will discuss the properties under finite
magnetic fields $\bm{h}\neq{\bf 0}$. It turns out that one needs
$O(\varepsilon)$ terms in the Zeeman energy.

###  III.1 Two-sublattice mean-field approximation for zero magnetic field

Throughout this paper we assume two-sublattice orders, which are indeed
observed experimentally in the double-perovskite compounds. This is because
our primary purposes in this paper are to demonstrate the electronic mechanism
for the induced ferro quadrupole moments under antiferro quadrupole orders
(Sec. III.1) and the finite magnetic moments in $J_{\rm eff}=3/2$ systems
(Sec. III.2). In recent analysis, Ref. Svoboda  _et al._ , 2021 have carried
out similar calculations with a four-sublattice unit cell. However, there is
no comment about the finite ferro quadrupole moments in the antiferroic
quadrupole phase, although the reason is unclear. For our purposes, it is
sufficient to use the two-sublattice unit cell. In some parameter sets used in
the following sections, the four-sublattice antiferro magnetic ordersSvoboda
_et al._ (2021), which correspond to double-$\bm{q}$ magnetic orders with
induced a single-$\bm{q}$ quadrupole moments, are more stabilized than FM110
states (see below), when one applies the 4-sublattice approximation. In this
sense, our results within the two-sublattice orders are regarded as simple
extension with the excited states from the $J_{\rm eff}=3/2$ model in Ref.
Chen  _et al._ , 2010. Nevertheless, the comparison between the results within
the two-sublattice approximation and the experiments for the FM110 orders in
Sec. III.2 possesses important information since regardless of which types of
the microscopic interactions favor the experimentally observed FM110 phase,
qualitative properties in the FM110 phase are unchanged. The detail analysis
about the four-sublattice orders is one of our future problems.

Let us take the two sublattice position in the unit cell as A: $(0,0,0)$ and
B: $(0,\frac{1}{2},\frac{1}{2})$ and assume the ordered configurations are
uniform on the $xy$ planes: $(x,y,n)$ and $(x+\frac{1}{2},y,n+\frac{1}{2})$,
where the lattice constant is set to unity and $n$ is an integer and $x$ and
$y$ represent the site positions on the $z=n$ plane. This choice corresponds
to the domain with its ordering wavevector ${\bm{k}}^{*}_{3}\equiv(0,0,2\pi)$.
The expressions for the mean-field Hamiltonian are trivially obtained and
listed in Appendix C. The inclusion of the CEF excited doublet in the model
Hamiltonian causes only small effects on the overall feature of the phase
diagram and the magnitude of the primary order parameters. Thus, the
qualitative results here are almost identical to those in Ref. Chen  _et al._
, 2010 except for the ferroic quadrupole moments.

Figure 2 shows the temperature dependence of the ferro and antiferro
quadrupole moments defined as

$\displaystyle{\bm{Q}}^{\rm F}\equiv\langle{\bm{Q}}\rangle_{\rm
A}+\langle{\bm{Q}}\rangle_{\rm B},\quad{\bm{Q}}^{\rm
AF}\equiv\langle{\bm{Q}}\rangle_{\rm A}-\langle{\bm{Q}}\rangle_{\rm B}.$ (28)

Here $\langle\cdot\rangle_{\rm A(B)}$ indicates the expectation value
calculated by the A(B)-site local mean-field Hamiltonian. The data for
$\lambda/J=10^{4}\gg 1$ corresponds to those shown in Fig. 8 in Ref. Chen  _et
al._ , 2010. For these parameter sets, the system undergoes a phase transition
into a pure quadrupolar phase with the primary antiferroic and induced ferroic
ones (AFQ+$f$q) at $T/J\simeq 0.9$ and another transition at $T/J\simeq 0.3$
into a magnetic phase (FM110). The order parameter configuration of the FM110
phase is schematically shown in Fig. 2: the $x$ and $y$ components of the
ferro/antiferro dipole ${\bm{S}}^{\rm F,AF}$ and ferro/antiferro octupole
${\bm{T}}_{\alpha,\beta}^{\rm F,AF}$ moments take finite values. The
definition of these multipole moments are given in Appendix B and we use the
definition of “F/AF” similarly to Eq. (28). Note that the realistic value of
$\lambda$ in Fig. 2 is e.g., $\lambda/J=100$ baring the magnetic transition
temperature $\sim$20 K in mind. One can see a finite $Q_{u}^{\rm
F}\lesssim-0.1$ inside the AFQ+$f$q phase for $0.3\lesssim T/J\lesssim 0.9$.
It is quite natural that as $\lambda$ decreases and thus the energy of the
excited doublet $\epsilon_{1/2}$ lowers, the ferro moment $|Q_{u}^{\rm F}|$
increases. It is interesting that the effect of the excited state is not
negligible even if the excited state energy is 100 times larger than the
transition temperature into the quadrupole order. The presence of the excited
doublet also affects the transition temperatures and the two transition
temperatures slightly increase upon decreasing $\lambda$ but these are minor
points.

Figure 2: Temperature dependence of the quadrupole moments $Q_{u}^{\rm F}$ and
$Q_{v}^{\rm AF}$ for $J^{\prime}/J=0.3,\ V/J=0.3$, and $\lambda/J=10^{4},\
500,\ 100,\ 50$, and $20$. $Q_{v}^{\rm F}=Q_{u}^{\rm AF}=0$ and are not shown.

Figure 3 shows $T$-$\lambda$ phase diagram for $V/J=0.3$ and
$J^{\prime}/J=0.1$–$0.5$. For these parameter sets, the AFQ+$f$q phase
robustly appears from the normal (paramagnetic) state via the second-order
transition. For $J^{\prime}/J=0.3$ and $0.5$, the ground state is the FM110,
whose configuration is shown in the inset of Fig. 2. For $J^{\prime}/J=0.1$,
the ground states change with varying $\lambda$. When $\lambda$ is large, the
planer antiferromagnetic (AFM) order takes place, where $S^{\rm AF}_{y}/S^{\rm
AF}_{x}=T_{\alpha,y}^{\rm AF}/T_{\alpha,x}^{\rm AF}=-T_{\beta,y}^{\rm
AF}/T_{\beta,x}^{\rm AF}$ and the other multipole moments are zero. As noted
in Ref. Chen  _et al._ , 2010, owing to an accidental degeneracy there is no
anisotropy and $S^{\rm AF}_{x,y}=0$ at $T=0$, while the $x$ and $y$ components
of ${\bm{T}}^{\rm AF}_{\alpha,\beta}$ are finite. This feature continues even
for finite $\lambda$. This is because the ground-state wavefunction at each
sublattice (A and B) in the AFM phase for finite $\lambda$ is exactly the same
as that at $\lambda=\infty$, which is the eigenstate of the quadrupole moment
$\langle\bm{Q}\rangle_{\rm A,B}=(-1+\delta,0)$. As is evident from Eq. (11),
the matrix $u$ does not have finite matrix elements between the excited states
and the fifth and sixth states, while $v$ in Eq. (12) does. In the language of
${\bm{S}}_{\rho}$, the ground state is characterized by
$\langle{\bm{S}}_{xy}\rangle_{\rm A,B}={\bf 0}$, $\langle
S_{yz,zx}^{z}\rangle_{\rm A,B}=0$, $\langle S_{yz}^{\mu}\rangle_{\rm
A}=-\langle S_{yz}^{\mu}\rangle_{\rm B}=-\langle S_{zx}^{\mu}\rangle_{\rm
A}=\langle S_{zx}^{\mu}\rangle_{\rm B}$ with $\mu=x$ or $y$. They ensure that
the state with $\langle\bm{Q}\rangle_{\rm A,B}=(-1+\delta,0)$ remains to be
decoupled from the excited states since there is only a term proportional to
$S_{i,yz}^{\mu}-S_{i,zx}^{\mu}(\mu=x,y)$ in the magnetic sector of the mean-
field Hamiltonian at the site $i$. See Appendix C. This has indeed no matrix
element between the states with $u=-1+\delta$ and the excited states as in
$u$. See their expressions in Eqs. (62), (63), (65), and (66).

In contrast, for small $\lambda$, the low-$T$ phase is another
antiferromagnetic phase (AFM*) as schematically shown in the inset of Fig.
3(a). The order parameters for the AFM* phase are similar to those in the AFM
phase, but there are finite $S_{x,y}^{\rm AF}$ moments even at $T=0$. To
realize the AFM* phase, the spin-orbit coupling $\lambda$ is tuned to be
$\lambda\lesssim 30J$, where a finite quadrupole moments $v$ appears and thus
the symmetry is lowered. This also causes a finite mixing between the states
with $u=-1+\delta$ and the excited states. Thus, the ground-state wavefunction
in the AFM state is no longer eigen state of the mean-field Hamiltonian. Then,
a different state becomes the ground state with a finite ${\bm{S}}^{\rm AF}$.
Indeed, such transition has been studied in the non-Kramers $\Gamma_{3}$
system in our previous studyHattori  _et al._ (2022), where a kind of
topological protection is important. However, the value of the spin-orbit
coupling is not realistic, and thus, we do not analyze it in detail here.

Figure 3: $T$-$\lambda$ phase diagram for $V/J=0.3$, and (a)
$J^{\prime}/J=0.1$, (b) $0.3$, and (c) 0.5. Inset in (a): the
antiferromagnetic configurations ${\bm{S}}^{\rm AF}$, ${\bm{T}}_{\alpha}^{\rm
AF}$, and ${\bm{T}}_{\beta}^{\rm AF}$ in the AFM* phase.

Let us now examine the $T$-$V$ phase diagram for the realistic value of
$\lambda\simeq 100J$. There are three phases: the AFM, the AFQ+$f$q, and the
FM110 phases in Fig. 4. For visualizing the magnitude of the $Q^{\rm F}_{u}$
in the ordered phases, $\sqrt{-Q^{\rm F}_{u}}$ is also depicted in Fig. 4 as
colormap. The reason for using the “square-root” is just due to the technical
one to show the finiteness in the AFQ+$f$q phase, where the $|Q_{u}^{\rm F}|$
is much smaller than that in the other ordered phases. We have checked that
even for $\lambda/J=10000$ ignoring the excited states, the phase diagram is
semiquantitatively the same as that shown in Fig. 4. The difference is that
the $Q_{u}^{\rm F}\simeq 0$ in the AFQ+$f$q for $\lambda/J=10000$, which
corresponds to the results in Ref. Chen  _et al._ , 2010.

Figure 4: $T$-$V$ phase diagram for $J^{\prime}/J=0.2$. The lines represent
the phase boundaries. The color map shows $\sqrt{-Q_{u}^{\rm F}}$, where
$Q_{u}^{\rm F}\leq 0$ for any $V$ and $T$. The square-root of $-Q_{u}^{\rm F}$
is used just for visualization of the small value of $Q_{u}^{\rm F}$ in the
AFQ phase.

###  III.2 Finite magnetic fields

Under the magnetic field $\bm{h}=(h_{x},h_{y},h_{z})$, which includes the Bohr
magneton $\mu_{\rm B}$ and is related to the real magnetic field $\bm{H}$ as
$\bm{h}=\mu_{\rm B}\bm{H}$, we introduce the Zeeman coupling with the
electronic $g$ factor with $g=2$ as

$\displaystyle H_{\rm
Z}=-\sum_{\sigma,\sigma^{\prime},m,m^{\prime}}{\bm{h}}\cdot({\bm{M}})_{\sigma\sigma^{\prime},mm^{\prime}}d_{m\sigma}^{\dagger}d_{m^{\prime}\sigma^{\prime}}.$
(29)

Here, the magnetic moment operator is defined as
$({\bm{M}})_{\sigma\sigma^{\prime},mm^{\prime}}\equiv[(\mathbbm{1})_{\sigma\sigma^{\prime}}({\bm{L}})_{mm^{\prime}}+2\bm{s}_{\sigma\sigma^{\prime}}(\mathbbm{1})_{mm^{\prime}}]$.
In actual calculations for the six states manifold, $\bm{M}$ reduces to the
matrix form listed in Eqs. (56)–(58) and Eqs. (68)–(70).

To determine which ordering wavevectors are stabilized under the magnetic
fields within the two-sublattice mean-field approximation, we repeat the
calculations with different ordering wavevector set up:
${\bm{k}}^{*}_{1}=(2\pi,0,0)$, ${\bm{k}}^{*}_{2}=(0,2\pi,0)$, and
${\bm{k}}^{*}_{3}=(0,0,2\pi)$ and compare their free energies. We note that
the diagonal matrix elements of ${\bm{M}}$ in the bases of the eigenstates of
the Hamiltonian (1) are $O(\varepsilon)$ in the CEF ground-state quartet. For
example, the diagonal elements of $M_{z}$ is
diag$(1,-1,\frac{4}{5}\varepsilon+\frac{1}{50}\varepsilon^{2},-\frac{4}{5}\varepsilon-\frac{1}{50}\varepsilon^{2},0,0)+O(\varepsilon^{3})$.
Thus, the half of the quartet state with $u\simeq 1$ has a finite dipole
moment $\sim\frac{4}{5}\varepsilon\mu_{\rm B}\simeq 0.4\mu_{\rm B}$, which has
not been recognized well so far.

Before discussing the numerical data, we should note the followings. Using the
$O(\varepsilon)$ terms in the $\bm{M}$ means that we treat the CEF ground
state quartet as effective $t_{2g}$ state with hybridizing to the excited
quartet. This treatment is valid if the interaction (27) are not modified
significantly when the processes related to the $e_{g}$ electrons are also
taken account. In this study, we do not estimate such processes and we keep
the interaction form derived within the $t_{2g}$ states as a starting
approximation. The more sophisticated analysis remains as one of our future
problems.

Figure 5 shows the magnetization
$M_{h}\equiv\langle\bm{M}\rangle\cdot\bm{h}/h$ and $|\langle\bm{M}\rangle|$
for $J^{\prime}/J=0.3$ and $V/J=0.27$ as a function of $h=|\bm{h}|$ for three
high-symmetric directions $\bm{h}\parallel[001],\ [110]$, and $[111]$. At
$h=0$, the ground state is the FM110 phase and the intermediate-temperature
phase is the AFQ+$f$q phase for this parameter set. As for the scale of the
magnetic field, $h/J=0.1$ corresponds to $\sim 10$ T if we set the magnetic
transition temperature $T_{m}\sim 0.34J\sim 20$ K. There are three equivalent
domains for the FM110-type phases: FM110, FM101, and FM011. As can be
trivially understood, these phases have its magnetization along [110], [101],
and [011] directions, respectively. A finite magnetic field selects some of
the three domains and the results are shown only for those with the lowest
free energy.

Owing to the difference in the field direction and that for the spontaneous
one for $\bm{h}\parallel[001]$ and [111], $M\neq M_{h}$, while $M=M_{h}$ for
$\bm{h}\parallel[110]$. One can notice that the increase in $M_{h}$ as
increasing $h$ is steepest for $\bm{h}\parallel[001]$. This is related to the
fact that the ferro quadrupole moment $Q_{u}^{\rm F}$ is large in the FM110
phase (Fig. 2) and the increase in $M_{h}=M_{z}$ is more favored than that for
$\bm{h}\parallel[111]$. In the AFQ+$f$q phase, the magnetic field also selects
the two of the domains: $\bm{k}_{1}$ or $\bm{k}_{2}$. However, the energy
difference between the domains with $\bm{k}_{1,2}$ and $\bm{k}_{3}$ is so tiny
and it is not physically important in actual situations, where there are many
aspects such as electron-lattice couplings.

Figure 5: Magnetization curve for $J^{\prime}/J=0.3$ and $V/J=0.27$. (a)
${\bm{h}}\parallel[001]$, (b) ${\bm{h}}\parallel[110]$, and (c)
${\bm{h}}\parallel[111]$. The magnetization $M_{h}$ is defined as
$M_{h}\equiv\langle\bm{M}\rangle\cdot\bm{h}/h$, with $h\equiv|\bm{h}|$ and
plotted by lines, while the absolute value $M\equiv|\langle\bm{M}\rangle|$ is
indicated by symbols. The domain with the lowest free energy is considered;
(a) ${\bm{k}}_{1}$ or ${\bm{k}}_{2}$ domains, (b) ${\bm{k}}_{3}$ domain, and
(c) ${\bm{k}}_{1}$, ${\bm{k}}_{2}$, or ${\bm{k}}_{3}$ domains. The
temperatures shown correspond to the FM110 phase ($T/J=0.1<T_{m}/J\simeq
0.34$), the AFQ+$f$q phase ($T/J=0.5$), and the normal phase ($T/J=1.0$).

##  IV Discussions

In the following subsections, we discuss two aspects. The first is the induced
ferro quadrupole moments under the antiferro quadrupole orders. The second is
the magnitude of the magnetic moment, which vanishes in the model of $J_{\rm
eff}=3/2$.

###  IV.1 Ferro quadrupole moments

Let us compare the results in Sec. III and the experimental data in the double
perovskite compound focusing on the ferro quadrupole moments observed in
Ba2MgReO6Hirai  _et al._ (2020). The relative magnitude of the antiferro
quadrupole $Q^{\rm AF}_{v}$ and the ferro quadrupole $Q^{\rm F}_{u}$ moments
can be indirectly estimated by observing O displacement
$(\varepsilon_{u},\varepsilon_{v})$ corresponding to $(u,v)$ in the ordered
phases. At $T=6$ K, i.e., inside the magnetically-ordered phase, Hirai et al.,
reported that there is about 0.4 % elongation of the oxygen octahedron along
the $z(c)$ axis in average. Here “average” means that the analysis without
inplane $v$ type displacement. With the further analysis including this
inplane displacements, they estimated that the ratio $v/u\sim Q^{\rm
AF}_{v}/Q^{\rm F}_{u}\sim 4$. Here, we have assumed a linear relation between
the oxygen displacement and the $d$-electron quadrupole moments, ignoring
anharmonic couplings. In the Supplimental Material in Ref. Hirai  _et al._ ,
2020, the data of average oxygen positions at $T=25$ K inside the quadrupolar
phase are also listed. One can estimates 0.18 % elongation of the oxygen
octahedron along the $z(c)$ axis in average. This naively leads to the ratio
$Q^{\rm AF}_{v}/Q^{\rm F}_{u}\sim 8$. Although there is ambiguity about the
quadrupole-displacement coupling, which can be anisotropic even in the first-
order in $Q_{u,v}$, the finite values of $Q_{u}^{\rm F}$ shown in Fig. 2 are
qualitatively consistent with this. Indeed, when setting $J\sim 50$ K and
$\lambda$$/J$$=50$–$100$, the AFQ order appears at $T_{q}\sim 40$ K and the
magnetic one at $T_{m}\sim 25$ K in Fig. 2. In actual situation, one should
take into accoount couplings between the electrons and the oxygen displacement
as discussed recentlyIwahara and Chibotaru (2022); Iwahara  _et al._ (2018).
Nevertheless, this result demonstrates that the contribution of the excited
states on the induced ferro quadrupole moments in the AFQ state is a
noticeable in the double perovskite materials.

###  IV.2 Effective magnetic moment

We now discuss the magnetic properties of this system. It is well known that
the magnetic moment $\bm{M}$ in the effective $J_{\rm eff}=3/2$ theory, i.e.,
$\varepsilon\to 0$ limit, vanishes. This is one of the weak points in the
theory and some discussions about the impact of the orbital orders on the
magnetic susceptibility have been done recentlySvoboda  _et al._ (2021). Let
us estimate the effective moment $\mu_{\rm eff}$ defined by the high-
temperature asymptotic form of the magnetic susceptibility as

$\displaystyle\mu_{\rm B}M_{z}$ $\displaystyle\simeq$ $\displaystyle\mu_{\rm
B}\frac{\frac{4\varepsilon}{5}(e^{\frac{4}{5}\varepsilon
h/T}-e^{-\frac{4}{5}\varepsilon h/T})}{e^{\frac{4}{5}\varepsilon
h/T}+e^{-\frac{4}{5}\varepsilon
h/T}+2}\simeq\frac{(\frac{\sqrt{3}}{\sqrt{2}}\frac{4}{5}\varepsilon\mu_{\rm
B})^{2}}{3T}\frac{h}{\mu_{\rm B}},\ \ \ \ \ \ $ (30)

where we have assumed $\bm{h}\parallel[001]$. Thus,

$\displaystyle\mu_{\rm eff}=\frac{4\sqrt{3}}{5\sqrt{2}}\varepsilon\mu_{\rm
B}=\frac{4\sqrt{2}}{\sqrt{3}}\left(\frac{\frac{3}{2}\lambda}{5B_{4}}\right)\mu_{\rm
B}.$ (31)

As is evident from the factor $\frac{3}{2}\lambda/(5B_{4})$, this arises from
the hybridization between the ground and excited quartets. Remember that the
induced ferro quadrupole moments discussed in the previous subsection arises
from the offdiagonal elements of the quadrupole operators between the ground
state quartet and the first excited doublet. Substituting the values for
Rb2TaCl6 with $D\simeq 5B_{4}\simeq 3.2$ eV and $\lambda\simeq 0.27$ eV
leading to $\varepsilon\simeq 0.42$, the effective moment is estimated as
$\mu_{\rm B}\simeq 0.41\mu_{\rm B}$, which is similar to the observed one
$0.27\mu_{\rm B}$Ishikawa  _et al._ (2019). When using the typical values for
the double-perovskites, $3\lambda/2\sim 0.5$ eV Ahn  _et al._ (2017);
Paramekanti  _et al._ (2018) and $D\sim 5$ eVOh  _et al._ (2018), one finds
$\mu_{\rm eff}\simeq 0.49\mu_{\rm B}$. The present estimation is not far from
the observed one $\mu_{\rm eff}\simeq 0.68\mu_{\rm B}$Hirai and Hiroi (2019)
in Ba2MgReO6. The remaining discrepancy might be improved by taking into
account the effects of ligand ionsAhn  _et al._ (2017) and/or the vibronic
degrees of freedomIwahara and Chibotaru (2022). The fluctuation effects
neglected in the present analysis also affect the magnitude of the magnetic
moment quantitatively. Analyses including these more sophisticated aspects are
one of our future problems.

Compared to the estimation of $\mu_{\rm eff}$, the results of the ordered
magnetic moment shown in Fig. 5 is not satisfactory in the present mean-field
analysis. As discussed in Fig. 5, which roughly corresponds to the parameter
set for Ba2MgReO6 with setting $J\sim 50$ K, $M_{h}$ and $M$ are $\sim 0.1$.
This is only $\sim 20$–30 % of the observed moment in Ba2MgReO6Hirai and Hiroi
(2019). As shown in Fig. 2, $Q_{u}^{\rm F}<0$ and this corresponds to larger
occupation in the $yz$ and $zx$ orbitals at $h=0$. In terms of the four states
in the ground state quartet, the energy of the $\Downarrow$ states
($J^{z}_{\rm eff}=\pm 3/2$) is lower than that for the $\Uparrow$ states
($J^{z}_{\rm eff}=\pm 1/2$). See Eqs. (34)–(37) and the matrix $u$ in Eq.
(11). The point is that these $\Downarrow$ states have no moment even with the
$O(\varepsilon^{2})$ correction. Thus, the ordered moment is smaller compared
with the effective moment even without the additional $\sqrt{3}$ factor in the
definition of $\mu_{\rm eff}$: $0.49/\sqrt{3}=0.28>0.1$. Recently, Zhang et
al., have proposed that the quadrupole configurations in the AFQ+$f$q phase
modifies the spin-spin exchange interactions and the Dzyaloshinskii-Moriya
type interaction induced in the AFQ+$f$q phase stabilizes the magnetic order
of FM110 and obtained a reasonable value of the ferro moments,Zhang  _et al._
(2023) although the relation between the CEF excited quartet states are not
clear.

In Ref. Ahn  _et al._ , 2017, the orbital angular momentum renormalization was
discussed. Owing to the extended $d$ like orbital including the surrounding
$p$ orbitals at the oxygen sites, ${\bm{M}}=2{\bm{S}}+{\bm{L}}\to
2{\bm{S}}+\gamma{\bm{L}}$ with $\gamma=0.536$ for Ba2MgReO6. Although it is
unclear whether their analysis without the cubic anisotropy includes the
effect of the local excited states shown in this paper, let us here
qualitatively examine how this renormalization effect influences the results
in our model. By introducing $\gamma$ in Eq. (29), the value of $M$ and
$M_{h}$ change as decreasing $\gamma$ from $\gamma=1$. For the parameters
shown in Fig. 5, firstly $M$ and $M_{h}$ decrease and vanish at $\gamma\sim
0.8$ and then increase. At $\gamma=0.5$, $M_{h}\sim 0.25$ as shown in Fig. 6.
In actual situation in real materials, not only the Zeeman energy (29), but
also the interactions (27) must be modified. Nevertheless, a finite
$\varepsilon$ modifies the orbital character in the presence of the
renormalization owing to $\gamma$. As shown in Fig. 6, the orbital
contributing to the magnetic moment is $\sim\Uparrow$ state ($J^{z}_{\rm
eff}=\pm 1/2$) for $\gamma\sim 1$, while $\Downarrow$ components ($J_{\rm
eff}^{z}=\pm 3/2$) increase as lowering $\gamma$. For the ferro quadrupolar
phase observed in Rb2TaCl6 and Cs2TaCl6, $Q_{u}^{\rm F}>0$ and $\gamma\sim
0.8$ is estimatedIshikawa  _et al._ (2019). Since $Q_{u}^{\rm F}>0$,
$\Uparrow$ components (i.e., $J_{\rm eff}=\pm 1/2$) are energetically favored
both by the uniform distortions and magnetic fields.

For smaller $\gamma$, the correction owing to the excited state proportional
to $\varepsilon$ is minor in the present simplified analysis. However, for
large $\varepsilon$, $\mu_{\rm eff}\sim 0.5\mu_{\rm B}$ without the correction
owing to the ligand ions. It is expected that the combined corrections by the
excited states and the ligand ions influence the orbital character of the
magnetic moment in more complex ways than discussed here. The theory
describing the magnetic moment in the $J_{\rm eff}=3/2$ model requires more
elaborate treatments and the interpretation of the experimental results should
also be reconsidered. In this respect, it is highly important to clarify the
orbital profile, i.e., the ratio between $\Uparrow$ and $\Downarrow$, of the
magnetic moments as analyzed in, e.g., LiV2O4Shimizu  _et al._ (2012).

Figure 6: $\gamma$ dependence of $M,\ M_{h}$ for $T/J=0.1$ K, $h/J=0.05$ along
$\bm{h}\parallel[001]$, $J^{\prime}/J=0.3$ and $V/J=0.27$. The effective
moment $\mu_{\rm eff}$’s given by Eq. (31) for several values of $\varepsilon$
are also drawn.

##  V Summary

We have clarified that the ferroic components of $O_{20}$ are induced under
the antiferroic $O_{22}$ orders in the fcc lattice model. We have
microscopically demonstrated the mechanism of the induced ferroic moments.
This is in fact quite simple: just taking into account the CEF excited doublet
state. It is important that the excited state located at $\sim 5000$ K, $\sim
100$ times larger scale than the transition temperature $T_{q}\sim 30$ K, can
influence the order parameter via generating the local anisotropy. We have
also clarified that the magnetic dipole moment in the $J_{\rm eff}=3/2$ model
does not vanish in the realistic values of the spin-orbit coupling and
crystalline electric field, when the CEF excited quartet states are taken into
account. Since these aspects have not been seriously considered so far, we
believe our results shed a renewed light on the orbital orders in related
correlated systems.

## Acknowledgment

This work was supported by JSPS KAKENHI (Grant No. JP21H01031 and No.
JP23H04869).

##  Appendix A Wavefunctions

We list the lowest 6 local eigenstates for $\lambda>0$ and $B_{4}>0$. They are
mainly constructed by the $t_{2g}$ orbitals in the parameter regime relevent
to our considerations. In terms of $m\sigma$ bases
$\\{2\uparrow,1\uparrow,0\uparrow,-1\uparrow,-2\uparrow,2\downarrow,1\downarrow,0\downarrow,-1\downarrow,-2\downarrow\\}$,
the local eigenstates for the $\Gamma_{7}$ with the energy $\epsilon_{1/2}$
and the $\Gamma_{8}$ with $\epsilon_{3/2}$ are expressed as

$\displaystyle|\Gamma_{7}+\rangle:\ $
$\displaystyle\\{0,-2,0,0,0,-1,0,0,0,1\\},$ (32)
$\displaystyle|\Gamma_{7}-\rangle:\ $
$\displaystyle\left\\{-1,0,0,0,1,0,0,0,2,0\right\\},$ (33)
$\displaystyle|\Gamma_{8}+\Uparrow\rangle:\ $
$\displaystyle\left\\{0,\frac{\varepsilon}{4}+\frac{w}{4}+\frac{1}{2},0,0,0,-\frac{\varepsilon}{2}-\frac{w}{2},0,0,0,1\right\\},$
(34) $\displaystyle|\Gamma_{8}-\Uparrow\rangle:\ $
$\displaystyle\left\\{1,0,0,0,-\frac{\varepsilon}{2}-\frac{w}{2},0,0,0,\frac{\varepsilon}{4}+\frac{w}{4}+\frac{1}{2},0\right\\},$
(35) $\displaystyle|\Gamma_{8}-\Downarrow\rangle:\ $
$\displaystyle\left\\{0,0,0,-1-\frac{\varepsilon}{10}-\frac{w}{2},0,0,0,\frac{\sqrt{6}}{5}\varepsilon,0,0\right\\},$
(36) $\displaystyle|\Gamma_{8}+\Downarrow\rangle:\ $
$\displaystyle\left\\{0,0,\frac{\sqrt{6}}{5}\varepsilon,0,0,0,-1-\frac{\varepsilon}{10}-\frac{w}{2},0,0,0\right\\},$
(37)

where we have not normalized the eigenstates and
$w\equiv\sqrt{\varepsilon^{2}+4\varepsilon/5+4}$. We have labeled these states
by the diagonal matrix elements of $J_{z}$ and $u$ for $\varepsilon\to 0$. The
label $\pm$ represents the sign of the diagonal matrix element $J_{z}$ in Eq.
(9), while the quadrupole part $\Uparrow(\Downarrow)$ in the $\Gamma_{8}$
states indicates the positive (negative) diagonal matrix element $u$ in Eq.
(11).

##  Appendix B Multipole operators

We here list various multipole operators appearing in the Hamiltonian (27). In
the following, we omit the site index $i$ for simplicity.

$\displaystyle S_{xy}^{\mu}$
$\displaystyle=\frac{1}{3}\left[\frac{1}{2}\hat{\sigma}^{\mu}+\frac{1}{2}\hat{\sigma}^{\mu}\hat{u}\right],$
(38) $\displaystyle S_{yz}^{\mu}$
$\displaystyle=\frac{1}{3}\left[\frac{1}{2}\hat{\sigma}^{\mu}+\frac{1}{2}\hat{\sigma}^{\mu}\left(-\frac{1}{2}\hat{u}+\frac{\sqrt{3}}{2}\hat{v}\right)\right],$
(39) $\displaystyle S_{zx}^{\mu}$
$\displaystyle=\frac{1}{3}\left[\frac{1}{2}\hat{\sigma}^{\mu}+\frac{1}{2}\hat{\sigma}^{\mu}\left(-\frac{1}{2}\hat{u}-\frac{\sqrt{3}}{2}\hat{v}\right)\right].$
(40)

Here, $\hat{\sigma}^{\mu}/2$ represents the spin operator of the $t_{2g}$
orbital: $S_{xy}^{\mu}+S_{yz}^{\mu}+S_{zx}^{\mu}$, and
$\frac{1}{2}\hat{\sigma}^{\mu}\hat{u}$ and
$\frac{1}{2}\hat{\sigma}^{\mu}\hat{v}$ are the spin-orbital composite ones:
$2S_{xy}^{\mu}-S_{yz}^{\mu}-S_{zx}^{\mu}$ and
$\sqrt{3}(S_{yz}^{\mu}-S_{zx}^{\mu})$, respectively. In terms of the
irreducible representation of the Oh group, the operators belonging to
$\Gamma_{4}$ are

$\displaystyle{\bm{S}}\equiv\frac{1}{2}\\{\hat{\sigma}^{x},\hat{\sigma}^{y},\hat{\sigma}^{z}\\},$
(41)
$\displaystyle{\bm{T}}_{\alpha}\equiv\frac{1}{2}\left\\{\hat{\sigma}^{x}\left(-\frac{1}{2}\hat{u}+\frac{\sqrt{3}}{2}\hat{v}\right),\hat{\sigma}^{y}\left(-\frac{1}{2}\hat{u}-\frac{\sqrt{3}}{2}\hat{v}\right),\hat{\sigma}^{z}\hat{u}\right\\},$
(42)

and those belonging to $\Gamma_{5}$ are

$\displaystyle{\bm{T}}_{\beta}\equiv\frac{1}{2}\left\\{\hat{\sigma}^{x}\left(-\frac{\sqrt{3}}{2}\hat{u}-\frac{1}{2}\hat{v}\right),\hat{\sigma}^{y}\left(\frac{\sqrt{3}}{2}\hat{u}-\frac{1}{2}\hat{v}\right),\hat{\sigma}^{z}\hat{v}\right\\}.$
(43)

Using these irreps., one finds

$\displaystyle S_{xy}^{x}$
$\displaystyle=\frac{1}{3}\left[S^{x}-\frac{1}{2}T_{\alpha}^{x}-\frac{\sqrt{3}}{2}T_{\beta}^{x}\right],$
(44) $\displaystyle S_{xy}^{y}$
$\displaystyle=\frac{1}{3}\left[S^{y}-\frac{1}{2}T_{\alpha}^{y}+\frac{\sqrt{3}}{2}T_{\beta}^{y}\right],$
(45) $\displaystyle S_{xy}^{z}$
$\displaystyle=\frac{1}{3}(S^{z}+T_{\alpha}^{z}),$ (46) $\displaystyle
S_{yz}^{x}$ $\displaystyle=\frac{1}{3}(S^{x}+T_{\alpha}^{x}),$ (47)
$\displaystyle S_{yz}^{y}$
$\displaystyle=\frac{1}{3}\left[S^{y}-\frac{1}{2}T_{\alpha}^{y}-\frac{\sqrt{3}}{2}T_{\beta}^{y}\right],$
(48) $\displaystyle S_{yz}^{z}$
$\displaystyle=\frac{1}{3}\left[S^{z}-\frac{1}{2}T_{\alpha}^{z}+\frac{\sqrt{3}}{2}T_{\beta}^{z}\right],$
(49) $\displaystyle S_{zx}^{x}$
$\displaystyle=\frac{1}{3}\left[S^{x}-\frac{1}{2}T_{\alpha}^{x}+\frac{\sqrt{3}}{2}T_{\beta}^{x}\right],$
(50) $\displaystyle S_{zx}^{y}$
$\displaystyle=\frac{1}{3}(S^{y}+T_{\alpha}^{y}),$ (51) $\displaystyle
S_{zx}^{z}$
$\displaystyle=\frac{1}{3}\left[S^{z}-\frac{1}{2}T_{\alpha}^{z}-\frac{\sqrt{3}}{2}T_{\beta}^{z}\right].$
(52)

##  Appendix C Mean-field Hamiltonian

We show the detail expression of the two-site mean-field Hamiltonian used in
the maintext. The total mean-field Hamiltonian is $H_{\rm int}^{\rm mf}=H_{\rm
int}^{\rm A}+H_{\rm int}^{\rm B}+C$, where $H_{\rm int}^{\rm A(B)}$ is the
mean-field Hamiltonian at the A(B) sublattice and $C$ is a constant depending
on the order parameters. Denoting the expectation values of order parameters
$X$ at the A(B) sublattice as $\langle X\rangle_{\rm A(B)}$, we obtain $H_{\rm
int}^{\rm A}=\sum_{i\in{\rm A}}(H_{{\rm int},i}^{{\rm A},xy}+H_{{\rm
int},i}^{{\rm A},yz}+H_{{\rm int},i}^{{\rm A},zx})$, where

$\displaystyle H_{{\rm int},i}^{{\rm A},xy}=$ $\displaystyle z_{\rm
nn}\Big{\\{}\langle\bm{Q}\rangle_{\rm A}\cdot\Big{(}\tilde{g}_{\rm
iso}\mathbbm{1}+\tilde{g}_{\rm
ani}{\mathsf{K}}^{3}\Big{)}\bm{Q}_{i}-J^{\prime}\langle\bm{S}_{xy}\rangle_{\rm
A}\cdot\bm{S}_{i}$
$\displaystyle+\Big{[}(J+2J^{\prime})\langle\bm{S}_{xy}\rangle_{\rm
A}-J^{\prime}\langle\bm{S}\rangle_{\rm A}\Big{]}\cdot\bm{S}_{i,xy}\Big{\\}},$
(53) $\displaystyle H_{{\rm int},i}^{{\rm A},yz}=$ $\displaystyle z_{\rm
nn}\Big{\\{}\langle\bm{Q}\rangle_{\rm B}\cdot\Big{(}\tilde{g}_{\rm
iso}\mathbbm{1}+\tilde{g}_{\rm
ani}{\mathsf{K}}^{1}\Big{)}\bm{Q}_{i}-J^{\prime}\langle\bm{S}_{yz}\rangle_{\rm
B}\cdot\bm{S}_{i}$
$\displaystyle+\Big{[}(J+2J^{\prime})\langle\bm{S}_{yz}\rangle_{\rm
B}-J^{\prime}\langle\bm{S}\rangle_{\rm B}\Big{]}\cdot\bm{S}_{i,yz}\Big{\\}},$
(54)

where $i$ obviously belongs to the A sublattice sites and $z_{\rm nn}=4$ is
the number of the nearest neighbor sites on each of the $xy,yz$, and $zx$
plane. Similarly to $H_{\rm int}^{{\rm A},yz}$, $H_{\rm int}^{{\rm A},zx}$ can
be obtained by replacing $yz\to zx$ and ${\mathsf{K}}^{1}\to{\mathsf{K}}^{2}$
in Eq. (54). The mean-field Hamiltonian at the B sublattice $H_{\rm int}^{\rm
B}$ is also trivially obtained by replacing A$\leftrightarrow$B in Eqs. (53)
and (54). The constant part $C$ is necessary for calculating the free energy
and this is given as

$\displaystyle\frac{C}{z_{\rm nn}N}=$ $\displaystyle\sum_{\alpha={\rm
A,B}}\Bigg{[}-\langle\bm{Q}\rangle_{\alpha}\frac{\tilde{g}_{\rm
iso}\mathbbm{1}+\tilde{g}_{\rm
ani}{\mathsf{K}}^{3}}{2}\cdot\langle\bm{Q}\rangle_{\alpha}$
$\displaystyle+J^{\prime}\langle\bm{S}_{xy}\rangle_{\alpha}\cdot\langle\bm{S}\rangle_{\alpha}-\frac{J+2J^{\prime}}{2}|\langle\bm{S}_{xy}\rangle_{\alpha}|^{2}\Bigg{]}$
$\displaystyle+$ $\displaystyle\Bigg{[}-\langle\bm{Q}\rangle_{\rm
A}\Big{(}\tilde{g}_{\rm iso}\mathbbm{1}+\tilde{g}_{\rm
ani}{\mathsf{K}}^{1}\Big{)}\cdot\langle\bm{Q}\rangle_{\rm B}$
$\displaystyle+J^{\prime}\Big{(}\langle\bm{S}_{yz}\rangle_{\rm
A}\cdot\langle\bm{S}\rangle_{\rm B}+\langle\bm{S}\rangle_{\rm
A}\cdot\langle\bm{S}_{yz}\rangle_{\rm B}\Big{)}$
$\displaystyle-(J+2J^{\prime})\langle\bm{S}_{yz}\rangle_{\rm
A}\cdot\langle\bm{S}_{yz}\rangle_{\rm B}+(yz\to zx)\Bigg{]}.$ (55)

Note that $N$ is the number of the unit cell (= the number of the A
sublattice) and remember ${\mathsf{K}}^{1}\to{\mathsf{K}}^{2}$ when $yz\to
zx$.

##  Appendix D Matrix form of multipole operators

The explicit 6$\times$6 matrix forms of the spin and angular momentum
operators are listed below. First, the spin operators ${\bm{S}}$’s are given
with
$\alpha=\frac{3\sqrt{2}}{100}\varepsilon^{2}+O(\varepsilon^{3})\simeq\sqrt{2}\delta$
in Eq. (11) as,

$\displaystyle S^{x}$
$\displaystyle\simeq\frac{1}{6}\small\begin{bmatrix}0&1&0&\\!\\!\\!\\!-\sqrt{2}+\alpha&0&\\!\\!\\!\\!-\sqrt{6}+\sqrt{3}\alpha\\\
&0&\sqrt{2}-\alpha&0&\\!\\!\\!\\!-\sqrt{6}+\sqrt{3}\alpha&0\\\
&&0&\\!\\!\\!\\!-2+2\sqrt{2}\alpha&0&\\!\\!\\!\\!\sqrt{3}-\sqrt{6}\alpha\\\
&&&0&\\!\\!\\!\\!-\sqrt{3}+\sqrt{6}\alpha&0\\\ &&&&0&0\\\ &&&&&0\\\
\end{bmatrix},$ (56) $\displaystyle S^{y}$
$\displaystyle\simeq\frac{i}{6}\small\begin{bmatrix}0&1&0&\\!\\!\\!\\!-\sqrt{2}+\alpha&0&\\!\\!\\!\\!\sqrt{6}-\sqrt{3}\alpha\\\
&0&-\sqrt{2}+\alpha&0&\\!\\!\\!\\!-\sqrt{6}+\sqrt{3}\alpha&0\\\
&&0&\\!\\!\\!\\!-2+2\sqrt{2}\alpha&0&\\!\\!\\!\\!-\sqrt{3}+\sqrt{6}\alpha\\\
&&&0&\\!\\!\\!\\!-\sqrt{3}+\sqrt{6}\alpha&0\\\ &&&&0&0\\\ &&&&&0\\\
\end{bmatrix},$ (57) $\displaystyle S^{z}$
$\displaystyle\simeq\frac{1}{6}\small\begin{bmatrix}1&0&-2\sqrt{2}+2\alpha&0&0&0\\\
&-1&0&\\!\\!\\!\\!-2\sqrt{2}+2\alpha&0&0\\\ &&-1+\sqrt{2}\alpha&0&0&0\\\
&&&\\!\\!\\!\\!1-\sqrt{2}\alpha&0&0\\\ &&&&\\!\\!\\!\\!3-3\sqrt{2}\alpha&0\\\
&&&&&\\!\\!\\!\\!-3+3\sqrt{2}\alpha\\\ \end{bmatrix}.$ (58)

Here, we have not shown the left bottom part of the matrices since they are
all Hermitian matrices. For the $t_{2g}$ spin operators ${\bm{S}}_{\rho}$’s,

$\displaystyle S_{xy}^{x}$
$\displaystyle\simeq\frac{1}{6}\begin{bmatrix}0&1&0&-\sqrt{2}+\alpha&0&0\\\
&0&\sqrt{2}-\alpha&0&0&0\\\ &&0&-2+2\sqrt{2}\alpha&0&0\\\ &&&0&0&0\\\
&&&&0&0\\\ &&&&&0\\\ \end{bmatrix},$ (59) $\displaystyle S_{xy}^{y}$
$\displaystyle\simeq\frac{i}{6}\begin{bmatrix}0&1&0&-\sqrt{2}+\alpha&0&0\\\
&0&-\sqrt{2}+\alpha&0&0&0\\\ &&0&-2+2\sqrt{2}\alpha&0&0\\\ &&&0&0&0\\\
&&&&0&0\\\ &&&&&0\\\ \end{bmatrix},$ (60) $\displaystyle S_{xy}^{z}$
$\displaystyle\simeq\frac{1}{6}\begin{bmatrix}-1&0&-\sqrt{2}+\alpha&0&0&0\\\
&1&0&-\sqrt{2}+\alpha&0&0\\\ &&-2+2\sqrt{2}\alpha&0&0&0\\\
&&&2-2\sqrt{2}\alpha&0&0\\\ &&&&0&0\\\ &&&&&0\\\ \end{bmatrix},$ (61)
$\displaystyle S_{yz}^{x}$
$\displaystyle\simeq\frac{1}{6}\begin{bmatrix}0&-1&0&\frac{-1}{\sqrt{2}}+\frac{\alpha}{2}&0&\frac{-\sqrt{3}}{\sqrt{2}}+\frac{\sqrt{3}\alpha}{2}\\\
&0&\frac{1}{\sqrt{2}}-\frac{\alpha}{2}&0&\frac{-\sqrt{3}}{\sqrt{2}}+\frac{\sqrt{3}\alpha}{2}&0\\\
&&0&\frac{1}{2}-\frac{\alpha}{\sqrt{2}}&0&\frac{\sqrt{3}}{2}-\frac{\sqrt{3}\alpha}{\sqrt{2}}\\\
&&&0&\frac{-\sqrt{3}}{2}+\frac{\sqrt{3}\alpha}{\sqrt{2}}&0\\\
&&&&0&\frac{-3}{2}+\frac{3\alpha}{\sqrt{2}}\\\ &&&&&0\\\ \end{bmatrix},$ (62)
$\displaystyle S_{yz}^{y}$
$\displaystyle\simeq\frac{i}{6}\begin{bmatrix}0&1&0&\frac{1}{\sqrt{2}}-\frac{\alpha}{2}&0&\frac{\sqrt{3}}{\sqrt{2}}-\frac{\sqrt{3}\alpha}{2}\\\
&0&\frac{1}{\sqrt{2}}-\frac{\alpha}{2}&0&-\frac{\sqrt{3}}{\sqrt{2}}+\frac{\sqrt{3}\alpha}{2}&0\\\
&&0&\frac{-1}{2}+\frac{\alpha}{\sqrt{2}}&0&\frac{-\sqrt{3}}{2}+\frac{\sqrt{3}\alpha}{\sqrt{2}}\\\
&&&0&\frac{-\sqrt{3}}{2}+\frac{\sqrt{3}\alpha}{\sqrt{2}}&0\\\
&&&&0&\frac{3}{2}-\frac{3\alpha}{\sqrt{2}}\\\ &&&&&0\\\ \end{bmatrix},$ (63)
$\displaystyle S_{yz}^{z}$
$\displaystyle\simeq\frac{1}{6}\begin{bmatrix}1&0&\frac{-1}{\sqrt{2}}+\frac{\alpha}{2}&0&\frac{\sqrt{3}}{\sqrt{2}}-\frac{\sqrt{3}\alpha}{2}&0\\\
&-1&0&\frac{-1}{\sqrt{2}}+\frac{\alpha}{2}&0&\frac{-\sqrt{3}}{\sqrt{2}}+\frac{\sqrt{3}\alpha}{2}\\\
&&\frac{1}{2}-\frac{\alpha}{\sqrt{2}}&0&\frac{-\sqrt{3}}{2}+\frac{\sqrt{3}\alpha}{\sqrt{2}}&0\\\
&&&\frac{-1}{2}+\frac{\alpha}{\sqrt{2}}&0&\frac{-\sqrt{3}}{2}+\frac{\sqrt{3}\alpha}{\sqrt{2}}\\\
&&&&\frac{3}{2}-\frac{3\alpha}{\sqrt{2}}&0\\\
&&&&&\frac{-3}{2}+\frac{3\alpha}{\sqrt{2}}\\\ \end{bmatrix},$ (64)
$\displaystyle S_{zx}^{x}$
$\displaystyle\simeq\frac{1}{6}\begin{bmatrix}0&1&0&\frac{1}{\sqrt{2}}-\frac{\alpha}{2}&0&\frac{-\sqrt{3}}{\sqrt{2}}+\frac{\sqrt{3}\alpha}{2}\\\
&0&\frac{-1}{\sqrt{2}}+\frac{\alpha}{2}&0&\frac{-\sqrt{3}}{\sqrt{2}}+\frac{\sqrt{3}\alpha}{2}&0\\\
&&0&-\frac{1}{2}+\frac{\alpha}{\sqrt{2}}&0&\frac{\sqrt{3}}{2}-\frac{\sqrt{3}\alpha}{\sqrt{2}}\\\
&&&0&\frac{-\sqrt{3}}{2}+\frac{\sqrt{3}\alpha}{\sqrt{2}}&0\\\
&&&&0&\frac{3}{2}-\frac{3\alpha}{\sqrt{2}}\\\ &&&&&0\\\ \end{bmatrix},$ (65)
$\displaystyle S_{zx}^{y}$
$\displaystyle\simeq\frac{i}{6}\begin{bmatrix}0&-1&0&\frac{-1}{\sqrt{2}}+\frac{\alpha}{2}&0&\frac{\sqrt{3}}{\sqrt{2}}-\frac{\sqrt{3}\alpha}{2}\\\
&0&\frac{-1}{\sqrt{2}}+\frac{\alpha}{2}&0&-\frac{\sqrt{3}}{\sqrt{2}}+\frac{\sqrt{3}\alpha}{2}&0\\\
&&0&\frac{1}{2}-\frac{\alpha}{\sqrt{2}}&0&\frac{-\sqrt{3}}{2}+\frac{\sqrt{3}\alpha}{\sqrt{2}}\\\
&&&0&\frac{-\sqrt{3}}{2}+\frac{\sqrt{3}\alpha}{\sqrt{2}}&0\\\
&&&&0&\frac{-3}{2}+\frac{3\alpha}{\sqrt{2}}\\\ &&&&&0\\\ \end{bmatrix},$ (66)
$\displaystyle S_{zx}^{z}$
$\displaystyle\simeq\frac{1}{6}\begin{bmatrix}1&0&\frac{-1}{\sqrt{2}}+\frac{\alpha}{2}&0&-\frac{\sqrt{3}}{\sqrt{2}}+\frac{\sqrt{3}\alpha}{2}&0\\\
&-1&0&\frac{-1}{\sqrt{2}}+\frac{\alpha}{2}&0&\frac{\sqrt{3}}{\sqrt{2}}-\frac{\sqrt{3}\alpha}{2}\\\
&&\frac{1}{2}-\frac{\alpha}{\sqrt{2}}&0&\frac{\sqrt{3}}{2}-\frac{\sqrt{3}\alpha}{\sqrt{2}}&0\\\
&&&\frac{-1}{2}+\frac{\alpha}{\sqrt{2}}&0&\frac{\sqrt{3}}{2}-\frac{\sqrt{3}\alpha}{\sqrt{2}}\\\
&&&&\frac{3}{2}-\frac{3\alpha}{\sqrt{2}}&0\\\
&&&&&\frac{-3}{2}+\frac{3\alpha}{\sqrt{2}}\\\ \end{bmatrix}.$ (67)

Finally, we show the expression of the orbital angular momentum $\bm{L}$:

$\displaystyle L_{x}$
$\displaystyle\simeq\frac{1}{3}\small\begin{bmatrix}0&2&0&-\frac{1}{\sqrt{2}}+\frac{\alpha_{L}}{2}&0&-\frac{\sqrt{3}}{\sqrt{2}}+\frac{\sqrt{3}\alpha_{L}}{2}\\\
&0&\frac{1}{\sqrt{2}}-\frac{\alpha_{L}}{2}&0&-\frac{\sqrt{3}}{\sqrt{2}}+\frac{\sqrt{3}\alpha_{L}}{2}&0\\\
&&0&2-\eta_{L}&0&-\sqrt{3}+\zeta_{L}\\\ &&&0&\sqrt{3}-\zeta_{L}&0\\\
&&&&0&\kappa_{L}\\\ &&&&&0\\\ \end{bmatrix},$ (68)

$\displaystyle L_{y}$
$\displaystyle\simeq\frac{i}{3}\small\begin{bmatrix}0&2&0&-\frac{1}{\sqrt{2}}+\frac{\alpha_{L}}{2}&0&\frac{\sqrt{3}}{\sqrt{2}}-\frac{\sqrt{3}\alpha_{L}}{2}\\\
&0&-\frac{1}{\sqrt{2}}+\frac{\alpha_{L}}{2}&0&-\frac{\sqrt{3}}{\sqrt{2}}+\frac{\sqrt{3}\alpha_{L}}{2}&0\\\
&&0&2-\eta_{L}&0&\sqrt{3}-\zeta_{L}\\\ &&&0&\sqrt{3}-\zeta_{L}&0\\\
&&&&0&\kappa_{L}\\\ &&&&&0\\\ \end{bmatrix},$ (69) $\displaystyle L_{z}$
$\displaystyle\simeq\frac{1}{3}\small\begin{bmatrix}2&0&-\sqrt{2}+\alpha_{L}&0&0&0\\\
&-2&0&-\sqrt{2}+\alpha_{L}&0&0\\\ &&1+\beta_{L}&0&0&0\\\
&&&-1-\beta_{L}&0&0\\\ &&&&-3-\gamma_{L}&0\\\ &&&&&3+\gamma_{L}\\\
\end{bmatrix},$ (70)

where
$\kappa_{L}=\frac{9}{5}\varepsilon-\frac{9}{50}\varepsilon^{2}+O(\varepsilon^{3})$,
$\eta_{L}=\frac{3}{5}\varepsilon+\frac{3}{50}\varepsilon^{2}+O(\varepsilon^{3})$,
$\zeta_{L}=-\frac{3\sqrt{3}}{5}\varepsilon+\frac{3\sqrt{3}}{25}\varepsilon^{2}+O(\varepsilon^{3})$,
and
$\alpha_{L}=\frac{3\sqrt{2}}{5}\varepsilon-\frac{3}{50\sqrt{2}}\varepsilon^{2}+O(\varepsilon^{3})$,
$\beta_{L}=\frac{12}{5}\varepsilon-\frac{3}{10}\varepsilon^{2}-\frac{63}{250}\varepsilon^{3}+O(\varepsilon^{4})$,
and $\gamma_{L}=-\frac{9}{50}\varepsilon^{2}+O(\varepsilon^{3})$.

## References

  * Tokura and Nagaosa (2000) Y. Tokura and N. Nagaosa, Science 288, 462 (2000).
  * Fiebig  _et al._ (2016) M. Fiebig, T. Lottermoser, D. Meier, and M. Trassin, Nature Reviews Materials 1, 16046 (2016).
  * Kuramoto  _et al._ (2009) Y. Kuramoto, H. Kusunose, and A. Kiss, J. Phys. Soc. Jpn. 78, 072001 (2009).
  * Pesin and Balents (2010) D. Pesin and L. Balents, Nat. Phys. 6, 376 (2010).
  * Jackeli and Khaliullin (2009) G. Jackeli and G. Khaliullin, Phys. Rev. Lett. 102, 017205 (2009).
  * Nasu  _et al._ (2014) J. Nasu, M. Udagawa, and Y. Motome, Phys. Rev. Lett. 113, 197205 (2014).
  * Corboz  _et al._ (2012) P. Corboz, M. Lajkó, A. M. Läuchli, K. Penc, and F. Mila, Phys. Rev. X 2, 041013 (2012).
  * 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).
  * Nomoto  _et al._ (2016) T. Nomoto, K. Hattori, and H. Ikeda, Phys. Rev. B 94 (2016).
  * Hattori  _et al._ (2017) K. Hattori, T. Nomoto, T. Hotta, and H. Ikeda, J. Phys. Soc. Jpn. 86, 113702 (2017).
  * Watanabe  _et al._ (2013) H. Watanabe, T. Shirakawa, and S. Yunoki, Phys. Rev. Lett. 110, 027002 (2013).
  * Vasala and Karppinen (2015) S. Vasala and M. Karppinen, Prog. Solid State Chem. 43, 1 (2015).
  * Longo and Ward (1961) J. Longo and R. Ward, J. Am. Chem. Soc. 83, 2816 (1961).
  * Wiebe  _et al._ (2003) C. R. Wiebe, J. E. Greedan, P. P. Kyriakou, G. M. Luke, J. S. Gardner, A. Fukaya, I. M. Gat-Malureanu, P. L. Russo, A. T. Savici, and Y. J. Uemura, Phys. Rev. B 68, 134410 (2003).
  * Sleight  _et al._ (1962) A. W. Sleight, J. Longo, and R. Ward, Inorganic Chemistry 1, 245 (1962), https://doi.org/10.1021/ic50002a010 .
  * Chamberland and Levasseur (1979) B. L. Chamberland and G. Levasseur, Mater. Res. Bull. 14, 401 (1979).
  * Abakumov  _et al._ (1997) A. M. Abakumov, R. V. Shpanchenko, E. V. Antipov, O. I. Lebedev, and G. Van Tendeloo, J. Solid State Chem. 131, 305 (1997).
  * Wiebe  _et al._ (2002) C. R. Wiebe, J. E. Greedan, G. M. Luke, and J. S. Gardner, Phys. Rev. B 65 (2002).
  * Bramnik  _et al._ (2003) K. G. Bramnik, H. Ehrenberg, J. K. Dehn, and H. Fuess, Solid State Sci. 5, 235 (2003).
  * Yamamura  _et al._ (2006) K. Yamamura, M. Wakeshima, and Y. Hinatsu, J. Solid State Chem. 179, 605 (2006).
  * Marjerrison  _et al._ (2016) C. A. Marjerrison, C. M. Thompson, G. Sala, D. D. Maharaj, E. Kermarrec, Y. Cai, A. M. Hallas, M. N. Wilson, T. J. S. Munsie, G. E. Granroth, R. Flacau, J. E. Greedan, B. D. Gaulin, and G. M. Luke, Inorg. Chem. 55, 10701 (2016).
  * Hirai and Hiroi (2019) D. Hirai and Z. Hiroi, J. Phys. Soc. Jpn. 88, 064712 (2019).
  * Ishikawa  _et al._ (2019) H. Ishikawa, T. Takayama, R. K. Kremer, J. Nuss, R. Dinnebier, K. Kitagawa, K. Ishii, and H. Takagi, Phys. Rev. B 100, 045142 (2019).
  * Paramekanti  _et al._ (2018) A. Paramekanti, D. J. Singh, B. Yuan, D. Casa, A. Said, Y.-J. Kim, and A. D. Christianson, Phys. Rev. B 97, 235119 (2018).
  * Oh  _et al._ (2018) J. H. Oh, J. H. Kim, J. H. Jeong, and S. H. Chang, Curr. Appl. Phys. 18, 1225 (2018).
  * Paramekanti  _et al._ (2020) A. Paramekanti, D. D. Maharaj, and B. D. Gaulin, Phys. Rev. B 101, 054439 (2020).
  * Maharaj  _et al._ (2020) D. D. Maharaj, G. Sala, M. B. Stone, E. Kermarrec, C. Ritter, F. Fauth, C. A. Marjerrison, J. E. Greedan, A. Paramekanti, and B. D. Gaulin, Phys. Rev. Lett. 124, 087206 (2020).
  * Ahn  _et al._ (2017) K.-H. Ahn, K. Pajskr, K.-W. Lee, and J. Kuneš, Phys. Rev. B 95, 064416 (2017).
  * Chen  _et al._ (2010) G. Chen, R. Pereira, and L. Balents, Phys. Rev. B 82, 174440 (2010).
  * Lu  _et al._ (2017) L. Lu, M. Song, W. Liu, A. P. Reyes, P. Kuhns, H. O. Lee, I. R. Fisher, and V. F. Mitrović, Nat. Commun. 8, 14407 (2017).
  * Hirai  _et al._ (2020) D. Hirai, H. Sagayama, S. Gao, H. Ohsumi, G. Chen, T.-h. Arima, and Z. Hiroi, Phys. Rev. Research 2, 022063 (2020).
  * Iwahara and Chibotaru (2022) N. Iwahara and L. F. Chibotaru, Phys. Rev. B 107, L220404 (2023).
  * Hattori and Tsunetsugu (2014) K. Hattori and H. Tsunetsugu, J. Phys. Soc. Jpn. 83, 034709 (2014).
  * Hattori and Tsunetsugu (2016) K. Hattori and H. Tsunetsugu, J. Phys. Soc. Jpn. 85, 094001 (2016).
  * Ishitobi and Hattori (2021) T. Ishitobi and K. Hattori, Phys. Rev. B 104, L241110 (2021).
  * Tsunetsugu  _et al._ (2021) H. Tsunetsugu, T. Ishitobi, and K. Hattori, J. Phys. Soc. Jpn. 90, 043701 (2021).
  * Hattori  _et al._ (2022) K. Hattori, T. Ishitobi, and H. Tsunetsugu, Phys. Rev. B 107, 205126 (2023).
  * Ishitobi and Hattori (2023) T. Ishitobi and K. Hattori, Phys. Rev. B 107, 104413 (2023).
  * Svoboda  _et al._ (2021) C. Svoboda, W. Zhang, M. Randeria, and N. Trivedi, Phys. Rev. B 104, 024437 (2021).
  * Iwahara  _et al._ (2018) N. Iwahara, V. Vieru, and L. F. Chibotaru, Phys. Rev. B 98, 075138 (2018).
  * Zhang  _et al._ (2023) X. Zhang, J. Zou, and G. Xu, (2023), arXiv:2303.04010 [cond-mat.str-el] .
  * Shimizu  _et al._ (2012) Y. Shimizu, H. Takeda, M. Tanaka, M. Itoh, S. Niitaka, and H. Takagi, Nat. Commun. 3, 981 (2012).