# Multipolar spin liquid in an exactly solvable model for
$j_{\mathrm{eff}}=\frac{3}{2}$ moments

Vanuildo S. de Carvalho  Instituto de Física, Universidade Federal de Goiás,
74.001-970, Goiânia-GO, Brazil  Hermann Freire  Instituto de Física,
Universidade Federal de Goiás, 74.001-970, Goiânia-GO, Brazil  Rodrigo G.
Pereira  International Institute of Physics and Departamento de Física Teórica
e Experimental, Universidade Federal do Rio Grande do Norte, Campus
Universitário, Lagoa Nova, Natal, RN, 59078-970, Brazil

###### Abstract

We study an exactly solvable model with bond-directional quadrupolar and
octupolar interactions between spin-orbital entangled
$j_{\mathrm{eff}}=\frac{3}{2}$ moments on the honeycomb lattice. We show that
this model features a multipolar spin liquid phase with gapless fermionic
excitations. In the presence of perturbations that break time-reversal and
rotation symmetries, we find Abelian and non-Abelian topological phases in
which the Chern number evaluates to $0$, $\pm 1$, and $\pm 2$. We also
investigate quantum phase transitions out of the multipolar spin liquid using
a parton mean-field approach and orbital wave theory. In the regime of strong
integrability-breaking interactions, the multipolar spin liquid gives way to
ferroquadrupolar-vortex and antiferro-octupolar ordered phases that harbor a
hidden spin-$\frac{1}{2}$ Kitaev spin liquid. Our work unveils mechanisms for
unusual multipolar orders and quantum spin liquids in Mott insulators with
strong spin-orbit coupling.

##  I Introduction

Quantum spin liquids (QSLs) refer to disordered ground states of interacting
spin systems which are characterized by the emergence of fractionalized
excitations, topological order, and long-range entanglement Balents (2010);
Broholm  _et al._ (2020). The quest for QSLs stands as one of the most
prolific topics of investigation in condensed matter physics, as demonstrated
by the growing list of QSL candidates and by notable advances in our
theoretical understanding of phases beyond the Landau paradigm Savary and
Balents (2017); Takagi  _et al._ (2019); Trebst and Hickey (2022).

A guiding principle in the search for new materials is rooted in the
engineering of bond-dependent interactions between
$j_{\mathrm{eff}}=\frac{1}{2}$ moments, as found in the exactly solvable model
proposed by Kitaev for the honeycomb lattice Kitaev (2006). This model
exhibits a genuine QSL ground state resulting from the fractionalization of
$j_{\mathrm{eff}}=\frac{1}{2}$ moments into Majorana fermions and a
$\mathbb{Z}_{2}$ gauge field. Kitaev interactions appear as the leading term
in effective models for spin-orbit-assisted Mott insulators with $4d^{5}$ and
$5d^{5}$ electronic configuration Jackeli and Khaliullin (2009). Over the last
decade, several experiments have indicated the proximity of Kitaev material
candidates, such as the iridates (Na1-xLix)2IrO3 Singh and Gegenwart (2010);
Singh  _et al._ (2012); Cao  _et al._ (2013) and H3LiIr2O6 Kitagawa  _et al._
(2018), to a QSL phase. Special attention has been devoted to the compound
$\alpha$-RuCl3 Plumb  _et al._ (2014); Kasahara  _et al._ (2018); Tanaka  _et
al._ (2022), for which a half-quantized thermal Hall conductance has been
reported Kasahara  _et al._ (2018).

In recent years, the search for new systems with QSL ground states has
expanded to Mott insulators with effective magnetic moments beyond the
$j_{\mathrm{eff}}=\frac{1}{2}$ picture Witczak-Krempa  _et al._ (2014);
Takayama  _et al._ (2021). In that regard, higher-spin Kitaev models were
microscopically derived for the honeycomb Mott insulators NiI2 ($d^{8}$
configuration and spin $S=1$) Stavropoulos  _et al._ (2019) and CrI3 ($d^{3}$
configuration and $S=\frac{3}{2}$) Xu  _et al._ (2020, 2018); Stavropoulos
_et al._ (2021).

In contrast to the original spin-$\frac{1}{2}$ Kitaev model, these
generalizations are not integrable. Their properties have been studied by
semiclassical analysis Baskaran  _et al._ (2008), exact diagonalization
Stavropoulos  _et al._ (2019); Koga  _et al._ (2018); Hickey  _et al._ (2020),
the density matrix renormalization group method Dong and Sheng (2020); Zhu
_et al._ (2020); Khait  _et al._ (2021), and parton mean-field theory Jin  _et
al._ (2022). Remarkably, all these methods point to the existence of a QSL
ground state for both $S=1$ and $S=\frac{3}{2}$. In addition, Mott insulators
with $4d^{1}$ and $5d^{1}$ configurations have also attracted attention
Witczak-Krempa  _et al._ (2014); Takayama  _et al._ (2021). The effective spin
Hamiltonians for these systems are written in terms of
$j_{\mathrm{eff}}=\frac{3}{2}$ moments Chen  _et al._ (2010) and contain
multipolar interactions which can give rise to exotic hidden multipolar orders
and even spin-orbital liquid phases Natori  _et al._ (2016, 2017); Romhányi
_et al._ (2017); Natori  _et al._ (2018); Yamada  _et al._ (2018); Ishikawa
_et al._ (2019); de Farias  _et al._ (2020); Chen and Wu (2021); Natori  _et
al._ (2023).

In this work, we propose an exactly solvable model for
$j_{\mathrm{eff}}=\frac{3}{2}$ moments located on the sites of a honeycomb
lattice which could potentially describe some properties of $4d^{1}$ and
$5d^{1}$ materials. This model involves a combination of nearest-neighbor
bond-directional quadrupole-quadrupole and octupole-octupole interactions. The
representation of the quadrupolar and octupolar operators in terms of Majorana
fermions maps the model to noninteracting fermions hopping on the background
of a static $\mathbb{Z}_{2}$ gauge field. Since this model exhibits a gapless
quantum disordered ground state and spin-orbital correlations associated with
the fractionalization of higher magnetic multipoles, it epitomizes all the
features of a _multipolar spin liquid_ (MSL) Takayama  _et al._ (2021); Li and
Chen (2017); Sibille  _et al._ (2020); Rayyan  _et al._ (2023).

Our model should be viewed as an integrable point in the parameter space for
the most general model describing $j_{\mathrm{eff}}=\frac{3}{2}$ spin-orbital
systems. In this sense, the existence of the MSL offers a different avenue of
research in the quest for materials with unusual QSL properties. Indeed, we
show that, under the action of rotational and time-reversal-symmetry-breaking
fields, our exactly solvable model realizes non-trivial Abelian and non-
Abelian topological phases with gapless chiral Majorana modes on the edge.
Furthermore, we study the robustness of the MSL state against the formation of
long-range multipolar orders driven by integrability-breaking interactions.
Interestingly, even in the strong-coupling regime the system exhibits a hidden
QSL state described by a spin-$\frac{1}{2}$ Kitaev model on top of either a
ferroquadrupolar-vortex (FQV) or an antiferro-octupolar (AO) ordered state.

The remainder of our paper is structured as follows. In Sec. II, we introduce
and provide the motivation for the $j_{\mathrm{eff}}=\frac{3}{2}$ model with
bond-directional quadrupolar and octupolar interactions. In Sec. III, we study
the properties of the MSL in the exactly solvable limit of the model. While
the fully symmetric MSL is gapless, we show that magnetic and strain fields
can give rise to Abelian and non-Abelian topological phases. Section IV
focuses on the stability of the MSL and the quantum phase transitions out of
this state, analyzed using parton mean-field theory and orbital wave theory.
Finally, we present our concluding remarks in Sec. V. We leave to Appendixes A
and B some technical details on the identification of the lower-energy gauge-
flux sector of the MSL and also on the parton mean-field approach employed
here.

##  II Model

We consider interacting local moments with effective total angular momentum
$j_{\mathrm{eff}}=\frac{3}{2}$, which appear in Mott insulators with $4d^{1}$
or $5d^{1}$ transition-metal ions Chen  _et al._ (2010); Takayama  _et al._
(2021); Khomskii and Streltsov (2021). The single electron in the open shell
occupies the triply degenerate $t_{2g}$ orbitals with effective orbital
angular momentum $l_{\rm eff}=1$, which in turn are split by the spin-orbit
coupling into a higher-energy $j_{\rm eff}=1/2$ doublet and a lower-energy
$j_{\rm eff}=3/2$ quartet. The latter is formed by two Kramers doublets, which
may remain degenerate at low temperatures because the strong spin-orbit
coupling suppresses the Jahn-Teller effect Streltsov and Khomskii (2020).

An effective spin-orbital Hamiltonian for $j_{\mathrm{eff}}=\frac{3}{2}$
moments was derived in Refs. Yamada  _et al._ (2018); Natori  _et al._ (2018)
considering the dominant exchange process on the honeycomb lattice. The
starting point is a multiorbital Hubbard-Kanamori model that combines the
effects of the crystal electric field, spin-orbit coupling, and electronic
correlations, assuming a single hopping path for 90∘ bonds formed by edge-
sharing octahedra. The resulting Hamiltonian contains quadrupole-quadrupole
and octupole-octupole interactions with coupling constants of the same order
as the coupling between dipole moments. More generally, one should take into
account subleading hopping paths known to be important in the generic model
for Kitaev materials Rau  _et al._ (2014). Furthermore, the electrostatic
interaction between orbitals with different charge distributions and the
virtual exchange of optical phonons in the cooperative Jahn-Teller effect can
also renormalize the quadrupole-quadrupole interactions Baker (1971); Gehring
and Gehring (1975).

Out of the most general Hamiltonian for $j_{\mathrm{eff}}=\frac{3}{2}$
moments, here we shall focus on the regime of parameter space where the
leading interactions involve quadrupole and octupole moments. Such a hierarchy
of interactions could in principle be realized in layered materials where the
relative strength of the couplings can be tuned by strain Voleti  _et al._
(2021); Yadav  _et al._ (2018).

Table 1:  Quadrupolar and octupolar operators acting on $j_{\textrm{eff}}=\frac{3}{2}$ states. The overlines in $\overline{J^{\alpha}J^{\beta}}$ and $\overline{J^{\alpha}J^{\beta}J^{\gamma}}$ stand for a symmetrization with respect to the indices of the dipolar operators $J^{\alpha}$, e.g., $\overline{J^{x}(J^{y})^{2}}\equiv J^{x}(J^{y})^{2}+J^{y}J^{x}J^{y}+(J^{y})^{2}J^{x}$. Here we omit the octupolar operators that transform in the $\Gamma_{4}$ representation (see, e.g., Ref. Santini  _et al._ (2009)).  Moment | Symmetry | Operators
---|---|---
Quadrupole | $\Gamma_{3}$ | $O_{3}\equiv O^{3z^{2}-r^{2}}=\frac{1}{3}[3(J^{z})^{2}-\textbf{J}^{2}]$
|  | $O_{2}\equiv O^{x^{2}-y^{2}}=\frac{1}{\sqrt{3}}[(J^{x})^{2}-(J^{y})^{2}]$
| $\Gamma_{5}$ | $O^{x}\equiv O^{yz}=\frac{1}{\sqrt{3}}\overline{J^{y}J^{z}}$
|  | $O^{y}\equiv O^{zx}=\frac{1}{\sqrt{3}}\overline{J^{x}J^{z}}$
|  | $O^{z}\equiv O^{xy}=\frac{1}{\sqrt{3}}\overline{J^{x}J^{y}}$
Octupole | $\Gamma_{2}$ | $T^{xyz}=\frac{2}{3\sqrt{3}}\overline{J^{x}J^{y}J^{z}}$
| $\Gamma_{5}$ | $T^{x}\equiv T^{x}_{\beta}=\frac{2}{3\sqrt{3}}[\overline{J^{x}(J^{y})^{2}}-\overline{(J^{z})^{2}J^{x}}]$
|  | $T^{y}\equiv T^{y}_{\beta}=\frac{2}{3\sqrt{3}}[\overline{J^{y}(J^{z})^{2}}-\overline{(J^{x})^{2}J^{y}}]$
|  | $T^{z}\equiv T^{z}_{\beta}=\frac{2}{3\sqrt{3}}[\overline{J^{z}(J^{x})^{2}}-\overline{(J^{y})^{2}J^{z}}]$

We can label the multipolar operators acting on $j_{\mathrm{eff}}=\frac{3}{2}$
states by irreducible representations of the octahedral group Santini  _et
al._ (2009). As shown in Table 1, the quadrupole moment has three components,
$O^{\gamma}$, with $\gamma\in\\{1,2,3\\}\equiv\\{x,y,z\\}$, which forms a
$\Gamma_{5}$ representation, and two components, $O_{2}$ and $O_{3}$, which
transform in the $\Gamma_{3}$ representation. From the latter, we define the
linear combinations

$\tilde{O}^{\gamma}=\cos\left(\frac{2\pi\gamma}{3}\right)O_{3}+\sin\left(\frac{2\pi\gamma}{3}\right)O_{2},$
(1)

so that a $C_{3}$ rotation about the $(111)$ axis perpendicular to the lattice
plane acts as $\gamma\mapsto\gamma+1$ for both $O^{\gamma}$ and
$\tilde{O}^{\gamma}$. We then consider bond-dependent quadrupole-quadrupole
interactions of the form

$H_{q}=\sum_{\gamma=x,y,z}\sum_{\langle
jl\rangle_{\gamma}}(K_{q}O_{j}^{\gamma}O_{l}^{\gamma}+K_{q}^{\prime}\tilde{O}_{j}^{\gamma}\tilde{O}_{l}^{\gamma}),$
(2)

where $\langle jl\rangle_{\gamma}$ refers to nearest-neighbor sites along a
$\gamma$ bond of the honeycomb lattice, as in the standard notation for the
Kitaev model Trebst and Hickey (2022). The first term in $H_{q}$ can be
generated by integrating out vibronic couplings to trigonal distortions, which
are significant in some cubic $d^{1}$ systems Iwahara  _et al._ (2018);
Khomskii and Streltsov (2021). This term was considered in the model for a
quadrupolar spin liquid defined in Ref. de Farias  _et al._ (2020). The model
with the pure $K_{q}$ interaction has a macroscopically degenerate ground-
state manifold de Farias  _et al._ (2020). Adding perturbations to this
special point in the phase diagram, we can access different phases with exotic
multipolar orders or spin-orbital liquid behavior. On the other hand, the
second term in Eq. (2) can be generated by coupling to tetragonal distortions
or by electrostatic interactions Chen  _et al._ (2010). This term is known to
favor quadrupolar order with a vortex pattern Wu (2008); Khaliullin  _et al._
(2021).

We now add octupole-octupole interactions to our model. To reduce the number
of parameters, we focus on interactions that involve the components
transforming in the $\Gamma_{2}$ and $\Gamma_{5}$ representations shown in
Table 1. This means that we neglect the octupolar operators, which, like the
dipole moment $\mathbf{J}$, transform in the $\Gamma_{4}$ representation
Santini  _et al._ (2009). As we shall see in Sec. III, this choice is
convenient for the purpose of identifying exactly solvable limits of the
general model. We then consider

$H_{o}=\sum_{\gamma=x,y,z}\sum_{\langle
jl\rangle_{\gamma}}(K_{o}T_{j}^{\gamma}T_{l}^{\gamma}+K_{o}^{\prime}T^{xyz}_{j}T^{xyz}_{l}).$
(3)

Note that the octupolar operators are odd under time reversal, but the
interactions are time-reversal invariant. A bond-dependent octupolar
interaction similar to the first term in Eq. (3) appears in phenomenological
models for the hidden multipolar order in the heavy-fermion compound
CexLa1-xB6 Kubo and Kuramoto (2003, 2004). The second term in Eq. (3) favors
an octupolar order without lattice distortions analogous to the one proposed
for double perovskites with $j_{\rm eff}=2$ local moments Paramekanti  _et
al._ (2020); Churchill and Kee (2022). Altogether, the Hamiltonian
$H=H_{q}+H_{o}$ captures the competition between different types of multipolar
order.

To analyze the multipolar interactions, it is useful to express the
$j_{\mathrm{eff}}=\frac{3}{2}$ states in terms of pseudospin and a pseudo-
orbital degrees of freedom Shiina  _et al._ (1997); Natori  _et al._ (2016,
2017); de Farias  _et al._ (2020). The local Hilbert space is spanned by the
eigenstates $|m_{J}\rangle$ of $J^{z}$, which obey
$J^{z}|m_{J}\rangle=m_{J}|m_{J}\rangle$, with
$m_{J}\in\\{\pm\frac{1}{2},\pm\frac{3}{2}\\}$. We define the mapping

$\displaystyle\Big{|}m_{J}=\frac{3}{2}\Big{\rangle}$
$\displaystyle=\Big{|}s^{z}=-\frac{1}{2},\tau^{z}=\frac{1}{2}\Big{\rangle},$
(4) $\displaystyle\Big{|}m_{J}=\frac{1}{2}\Big{\rangle}$
$\displaystyle=-\Big{|}s^{z}=\frac{1}{2},\tau^{z}=-\frac{1}{2}\Big{\rangle},$
(5) $\displaystyle\Big{|}m_{J}=-\frac{1}{2}\Big{\rangle}$
$\displaystyle=\Big{|}s^{z}=-\frac{1}{2},\tau^{z}=-\frac{1}{2}\Big{\rangle},$
(6) $\displaystyle\Big{|}m_{J}=-\frac{3}{2}\Big{\rangle}$
$\displaystyle=-\Big{|}s^{z}=\frac{1}{2},\tau^{z}=\frac{1}{2}\Big{\rangle},$
(7)

where $\mathbf{s}_{j}$ and $\bm{\tau}_{j}$ are, respectively, pseudospin and
pseudo-orbital operators, which obey the SU(2) algebra
$[s^{\alpha}_{j},s^{\beta}_{k}]=i\epsilon^{\alpha\beta\gamma}\delta_{jk}s^{\gamma}_{j}$
and
$[\tau^{\alpha}_{j},\tau^{\beta}_{k}]=i\epsilon^{\alpha\beta\gamma}\delta_{jk}\tau^{\gamma}_{j}$
and the relation $[s^{\alpha}_{j},\tau^{\beta}_{k}]=0$. In the
$\\{|s^{z},\tau^{z}\rangle\\}$ basis, the quadrupolar and octupolar operators
given in Table 1 become

$\displaystyle O^{\gamma}$ $\displaystyle=$
$\displaystyle-4s^{\gamma}\tau^{y},\qquad O_{3}=2\tau^{z},\qquad
O_{2}=2\tau^{x},$ $\displaystyle T^{\gamma}$ $\displaystyle=$
$\displaystyle-4s^{\gamma}\mathbf{v}^{\gamma}\cdot\bm{\tau},\qquad
T^{xyz}=2\tau^{y},$ (8)

where
$\mathbf{v}^{\gamma}=\cos(2\pi\gamma/3)\hat{\mathbf{x}}-\sin(2\pi\gamma/3)\hat{\mathbf{z}}$
are unit vectors in the $xz$ plane of pseudo-orbital space.

##  III Exactly solvable model

In this section, we analyze the Hamiltonian $H=H_{q}+H_{o}$ in the special
limit $K_{q}^{\prime}=K_{o}^{\prime}=0$, in which we obtain

$H_{s}=\sum_{\gamma=x,y,z}\sum_{\langle
jl\rangle_{\gamma}}(K_{q}O_{j}^{\gamma}O_{l}^{\gamma}+K_{o}T_{j}^{\gamma}T_{l}^{\gamma}).$
(9)

We show that this model can be solved exactly using a Majorana fermion
representation for the multipolar operators and exhibits a MSL ground state.

Figure 1: Schematic representation of fractionalized
$j_{\mathrm{eff}}=\frac{3}{2}$ moments on the honeycomb lattice. Here $x$,
$y$, and $z$ bonds are represented by red, green, and blue solid lines,
respectively. The two sublattices, A and B, are represented by up-pointing and
down-pointing triangles. The inset in the top right corner shows the matter
($\theta^{\gamma}$) and gauge ($\eta^{\gamma}$) Majorana fermions. Each
plaquette is associated with the conserved quantity $\hat{W}_{\small
p}=-e^{i\pi(J^{x}_{1}+J^{y}_{2}+J^{z}_{3}+J^{x}_{4}+J^{y}_{5}+J^{z}_{6})}$.

We can write the pseudospin and pseudo-orbital operators in terms of six
Majorana fermions as Natori  _et al._ (2016); de Farias  _et al._ (2020)

$\displaystyle s_{j}^{\gamma}$
$\displaystyle=-\frac{i}{4}\epsilon^{\alpha\beta\gamma}\eta_{j}^{\alpha}\eta_{j}^{\beta},$
(10) $\displaystyle\tau_{j}^{\gamma}$
$\displaystyle=-\frac{i}{4}\epsilon^{\alpha\beta\gamma}\theta_{j}^{\alpha}\theta_{j}^{\beta},$
(11)

where $\epsilon^{\alpha\beta\gamma}$ is the Levi-Civita tensor and
$\eta^{\alpha}$ and $\theta^{\alpha}$ are Majorana fermion operators which
obey the algebra
$\\{\eta^{\alpha}_{j},\eta^{\beta}_{k}\\}=2\delta_{jk}\delta^{\alpha\beta}$,
$\\{\theta^{\alpha}_{j},\theta^{\beta}_{k}\\}=2\delta_{jk}\delta^{\alpha\beta}$,
and $\\{\eta^{\alpha}_{j},\theta^{\beta}_{k}\\}=0$. This Majorana
representation enlarges the Hilbert space, introducing unphysical states. To
avoid that, we impose the local constraint

$D_{j}\equiv
i\eta_{j}^{x}\eta_{j}^{y}\eta_{j}^{z}\theta_{j}^{x}\theta_{j}^{y}\theta_{j}^{z}=1$
(12)

for every site $j$. This is equivalent to obtaining a physical state
$|\Psi\rangle_{\mathrm{phys}}$ from the Majorana wave function
$|\Psi_{0}\rangle$ by implementing the projection
$|\Psi\rangle_{\mathrm{phys}}={\Pi_{j}[\frac{1}{2}(1+D_{j})]|\Psi_{0}\rangle}$.
Using the local constraint and the expressions in Eq. (8), we can write the
multipolar operators as Majorana bilinears:

$\displaystyle O^{\gamma}$ $\displaystyle=$ $\displaystyle
i\eta^{\gamma}\theta^{y},\qquad O_{3}=-i\theta^{x}\theta^{y},\qquad
O_{2}=-i\theta^{y}\theta^{z},$ $\displaystyle T^{\gamma}$ $\displaystyle=$
$\displaystyle-i\eta^{\gamma}(\mathbf{v}^{\gamma}\cdot\bm{\theta}),\qquad
T^{xyz}=-i\theta^{z}\theta^{x}.$ (13)

Note that
$\mathbf{v}^{\gamma}\cdot\bm{\theta}=\cos(2\pi\gamma/3)\theta^{x}-\sin(2\pi\gamma/3)\theta^{z}$
only involves only $\theta^{x}$ and $\theta^{z}$.

Figure 2: Dispersion relation for Majorana fermion excitations in the exactly
solvable MSL model. Here we set $K_{o}=0.4K_{q}$. (a) For the unperturbed
model in Eq. (9), the system exhibits a gapless spectrum with Dirac nodes at
the $\mathbf{K}$ point. (b) In the presence of strain fields $\varepsilon_{2}$
and $\varepsilon_{3}$ and time-reversal-symmetry-breaking field $h$ [see Eq.
(17)], the excitation spectrum becomes gapped. Here we set
$\varepsilon_{2}=\varepsilon_{3}=0.4K_{q}$ and $h=0.1K_{q}$.

Using the Majorana representation, we cast the Hamiltonian in the form

$H_{s}=i\sum_{\gamma}\sum_{\langle jl\rangle_{\gamma}}\hat{u}_{\langle
jl\rangle_{\gamma}}[K_{q}\theta^{y}_{j}\theta^{y}_{l}+K_{o}(\mathbf{v}^{\gamma}\cdot\bm{\theta}_{j})(\mathbf{v}^{\gamma}\cdot\bm{\theta}_{l})],$
(14)

where $\hat{u}_{\langle
jl\rangle_{\gamma}}\equiv-i\eta^{\gamma}_{j}\eta^{\gamma}_{l}$ are
antisymmetric bond operators obeying $\hat{u}_{\langle
jl\rangle_{\gamma}}=-\hat{u}_{\langle lj\rangle_{\gamma}}$. These operators
commute with each other and define $\mathbb{Z}_{2}$ gauge fields since $H_{s}$
is invariant with respect to the local transformations
$\theta^{\gamma}_{j}\mapsto\Lambda_{j}\theta^{\gamma}_{j}$ and
$\hat{u}_{\langle jl\rangle_{\gamma}}\mapsto\Lambda_{j}\hat{u}_{\langle
jl\rangle_{\gamma}}\Lambda_{l}$, with $\Lambda_{j}=\pm 1$. Like what happens
for the dipolar Kitaev models Kitaev (2006); Baskaran  _et al._ (2008), the
Hamiltonian in Eq. (9) possesses a conserved operator $\hat{W}_{\small
p}\equiv-e^{i\pi(J^{x}_{1}+J^{y}_{2}+J^{z}_{3}+J^{x}_{4}+J^{y}_{5}+J^{z}_{6})}$
for each hexagonal plaquette $p$. In the Majorana representation, the
plaquette operators read $\hat{W}_{\small p}=\prod_{\small\langle
j,l\rangle_{\gamma}\in p}\hat{u}_{\langle j,l\rangle_{\gamma}}$ de Farias  _et
al._ (2020); Jin  _et al._ (2022). Due to the form of the Hamiltonian in Eq.
(14), we refer to $\eta^{\gamma}$ as the gauge fermions and to
$\theta^{\gamma}$ as the matter fermions. The fractionalization of the
multipolar operators into Majorana fermions is illustrated in Fig. 1.

We verified numerically that the ground state of $H_{s}$ lies in the sector
where the eigenvalues of all conserved plaquette operators are set to
$W_{p}=+1$ (see Appendix A). Fixing a gauge configuration with ${u}_{\langle
jl\rangle_{\gamma}}=1$ for sites $j$ in sublattice A and $l$ in sublattice B,
we obtain from Eq. (14) a translation-invariant quadratic Hamiltonian for the
matter fermions moving in the background of the static $\mathbb{Z}_{2}$ gauge
field. Note that the $\theta^{y}$ fermion decouples from $\theta^{x}$ and
$\theta^{z}$. The dynamics of $\theta^{y}$ is governed by the quadrupolar
coupling $K_{q}$, whereas the dynamics of $\theta^{x}$ and $\theta^{z}$
depends on the octupolar coupling $K_{o}$. We diagonalize the Hamiltonian
using the Fourier transform

$\theta^{\gamma}_{b}(\mathbf{R})=\sqrt{\frac{2}{N}}\sum_{\mathbf{k}\in\frac{1}{2}(\mathrm{BZ})}[e^{-i\mathbf{k}\cdot\mathbf{R}}\theta^{\gamma}_{b}(\mathbf{k})+e^{i\mathbf{k}\cdot\mathbf{R}}\theta^{\gamma\dagger}_{b}(\mathbf{k})],$
(15)

where $\mathbf{R}$ is the position of the unit cell, $b\in\\{A,B\\}$ refers to
the sublattice index, $N$ is the number of sites, and the momentum sum runs
over half of the Brillouin zone (BZ). The bands associated with the
$\theta^{y}$ fermions have the dispersion relation
$E_{s}(\mathbf{k})=\pm|K_{q}g_{y}(\mathbf{k})|$, where
$g_{y}(\mathbf{k})=1+e^{i\mathbf{k}\cdot\mathbf{n}_{1}}+e^{i\mathbf{k}\cdot\mathbf{n}_{2}}$,
with lattice vectors $\mathbf{n}_{1}=\frac{\sqrt{3}}{2}(1,\sqrt{3})$ and
$\mathbf{n}_{2}=\frac{\sqrt{3}}{2}(-1,\sqrt{3})$. For the bands associated
with $\theta^{x}$ and $\theta^{z}$, we obtain

$E_{s}(\mathbf{k})=\pm\frac{|K_{o}|}{\sqrt{2}}\sqrt{|g_{z}(\mathbf{k})|^{2}+|g_{x}(\mathbf{k})|^{2}+2|g_{xz}(\mathbf{k})|^{2}\pm\sqrt{[|g_{z}(\mathbf{k})|^{2}-|g_{x}(\mathbf{k})|^{2}]^{2}+4|g_{x}(\mathbf{k})g^{*}_{xz}(\mathbf{k})+g^{*}_{z}(\mathbf{k})g_{xz}(\mathbf{k})|^{2}}},$
(16)

where we define the functions
$g_{z}(\mathbf{k})={\frac{3}{4}(e^{i\mathbf{k}\cdot\mathbf{n}_{1}}+e^{i\mathbf{k}\cdot\mathbf{n}_{2}})}$,
$g_{x}(\mathbf{k})=1+{\frac{1}{4}(e^{i\mathbf{k}\cdot\mathbf{n}_{1}}+e^{i\mathbf{k}\cdot\mathbf{n}_{2}})}$,
and
$g_{xz}(\mathbf{k})={\frac{\sqrt{3}}{4}(e^{i\mathbf{k}\cdot\mathbf{n}_{1}}-e^{i\mathbf{k}\cdot\mathbf{n}_{2}})}$.

The band structure of the Majorana fermion excitations is shown in Fig. 2(a).
The spectrum exhibits gapless modes with Dirac cones at the $\mathbf{K}$ point
in both sectors associated with $\theta^{y}$ fermions and mixed $\theta^{x}$
and $\theta^{z}$ fermions. The latter sector also displays a completely flat
gapped band. This peculiar feature holds for only the model with spatially
isotropic couplings; the gapped band becomes dispersive as soon as we
introduce different octupolar couplings $K_{o}^{(\gamma)}$ for each bond
direction $\gamma$. Importantly, by generating a dispersion for $\theta^{x}$
and $\theta^{z}$ fermions, the octupolar interaction lifts the macroscopic
degeneracy of the pure quadrupolar model de Farias  _et al._ (2020).
Consequently, our exactly solvable model stabilizes a unique MSL ground state.

In addition to gapless Majorana fermions, the spectrum of $H_{s}$ contains
vortex excitations (visons) associated with changing the eigenvalues of the
plaquette operators to $W_{p}=-1$. These gapped excitations are created in
pairs by the action of local operators which anticommute with the bond
operators $\hat{u}_{\langle jl\rangle_{\gamma}}$. Remarkably, the quadrupole
operators $O_{2}$ and $O_{3}$ and the octupolar operator $T^{xyz}$ do not
excite visons since they commute with the gauge fermions [see Eq. (13)]. As a
result, $O_{2}$, $O_{3}$, and $T^{xyz}$ are the only on-site operators
featuring a gapless spectrum in the corresponding dynamical structure factors
probed by resonant inelastic x-ray scattering Natori  _et al._ (2017). By
contrast, the dynamical spin structure factor for the pure Kitaev model
exhibits a flux gap Knolle  _et al._ (2014).

In analogy with the Kitaev model, we can drive the MSL to topological phases
by adding perturbations that break symmetries and gap out the fermionic
spectrum. We can break the $\mathbb{Z}_{3}$ symmetry (rotation in real and
spin-orbital space) and time-reversal symmetry while still preserving the
integrability of the model if we consider the perturbations

$\displaystyle\delta H_{s}$
$\displaystyle=-\sum_{j}(\varepsilon_{2}O_{2,j}+\varepsilon_{3}O_{3,j}+hT^{xyz}_{j}).$
(17)

Here $\varepsilon_{2}$ and $\varepsilon_{3}$ describe tetragonal strain fields
and $h$ induces an octupole moment, allowed by symmetry in the presence of an
external magnetic field. According to Eq. (13), the total Hamiltonian remains
quadratic in the matter fermions, but the strain fields couple $\theta^{y}$ to
$\theta^{x}$ and $\theta^{z}$. As shown in Fig. 2(b), we find a fully gapped
spectrum for generic values of $\varepsilon_{2}$, $\varepsilon_{3}$, and $h$.
However, there exist critical values of these parameters for which the
Majorana gap closes and then reopens, which is suggestive of topological phase
transitions.

Figure 3: Topological phase diagram of the MSL as a function of the symmetry-
breaking fields $\varepsilon_{2}$ and $h$. Here we set $K_{o}=0.2K_{q}$ and
$\varepsilon_{2}=\varepsilon_{3}$. The purple, yellow, and blue regions
correspond to Chern numbers 0, 1, and 2, respectively.

To investigate the topological nature of the gapped phases, we evaluate the
Chern number

$\mathcal{N}_{\text{Ch}}\equiv\frac{1}{2\pi}\sum\limits_{n}\int_{\text{BZ}}d^{2}\mathbf{k}\,\mathcal{B}_{n}(\mathbf{k}),$
(18)

where
${\mathcal{B}}_{n}(\mathbf{k})=[\bm{\nabla\times}\bm{\mathcal{A}}_{n}(\mathbf{k})]_{z}$
is the Berry curvature associated with the $n$th energy band;
$\bm{\mathcal{A}}_{n}(\mathbf{k})=i\langle
u_{n}(\mathbf{k})|\bm{\nabla}_{\mathbf{k}}|u_{n}(\mathbf{k})\rangle$ is the
Berry connection, with $|u_{n}(\mathbf{k})\rangle$ being the eigenvectors of
the Hamiltonian; and the sum on the right-hand-side runs over the occupied
bands. We determine $\mathcal{N}_{\text{Ch}}$ numerically following the
approach devised in Ref. Fukui  _et al._ (2005). The dependence of
$\mathcal{N}_{\text{Ch}}$ on the symmetry-breaking fields is shown in Fig. 3.
Remarkably, the interplay between strain and magnetic fields realizes five of
Kitaev’s 16-fold-way phases Kitaev (2006); Zhang  _et al._ (2020);
Chulliparambil  _et al._ (2020): a trivial phase with
$\mathcal{N}_{\text{Ch}}=0$, two non-Abelian topological phases with
$\mathcal{N}_{\text{Ch}}=\pm 1$, and two Abelian topological phases with
$\mathcal{N}_{\text{Ch}}=\pm 2$. For $0<|\varepsilon_{2}|,|\varepsilon_{3}|$
and $|h|\ll|K_{q}|$, the system is in an Abelian phase with
$\mathcal{N}_{\text{Ch}}=\pm 2$, depending on the sign of $h$. We can reach
the non-Abelian phase by increasing the strength of the symmetry-breaking
fields, whereas the trivial phase appears only for strong strain fields.
Varying the ratio $K_{o}/K_{q}$, we observe that enhancing the octupolar
interaction favors the Abelian phase. We have verified that the Majorana gap
closes at topological phase transitions where the Chern number changes. The
momentum at which the gap closes depends on the interactions and is not
attached to the $\mathbf{K}$ point. A direct transition from
$\mathcal{N}_{\text{Ch}}=\pm 2$ to $\mathcal{N}_{\text{Ch}}=0$ is possible by
fine tuning so that the gap closes simultaneously at two different points in
the Brillouin zone.

As a signature of the topological phases, we look for chiral edge states in a
cylinder geometry with zigzag edges. The result matches the Chern number, in
accordance with the bulk-boundary correspondence. Figure 4 shows the spectrum
with two chiral edge states in the regime where $\mathcal{N}_{\rm Ch}=\pm 2$.
The chiral edge states can be probed by the thermal Hall conductance
$\kappa_{xy}$, which for low enough temperatures behaves as
$\kappa_{xy}=\frac{\mathcal{N}_{\text{Ch}}}{2}\frac{\pi k^{2}_{B}}{6\hbar}T$.

Figure 4: Edge states of the MSL for a honeycomb lattice with zigzag edges.
Here we use a cylinder with $L=50$ unit cells along the open direction. We fix
the parameters as $K_{o}=h=\varepsilon_{2}=\varepsilon_{3}=0.2K_{q}$ in the
regime of $\mathcal{N}_{\text{Ch}}=\pm 2$ (see Fig. 3).

##  IV Multipolar orders from integrability-breaking interactions

In this section, we address the effect of the interactions $K_{q}^{\prime}$
and $K_{o}^{\prime}$ in Eqs. (2) and (3) on the ground state and elementary
excitations of the spin-orbital system. We define

$H_{I}=\sum_{\gamma}\sum_{\langle
jl\rangle_{\gamma}}(K^{\prime}_{q}\tilde{O}_{j}^{\gamma}\tilde{O}_{l}^{\gamma}+K_{o}^{\prime}T^{xyz}_{j}T^{xyz}_{l}),$
(19)

so that $H=H_{q}+H_{o}=H_{s}+H_{I}$. Hereafter we assume that all coupling
constants are positive. In the Majorana representation in Eq. (13), $H_{I}$
gives rise to interaction terms which spoil the integrability of the model.
For this reason, we resort to two complementary analytical approximations
which provide a qualitative picture of the excitation spectrum, namely parton
mean-field theory and orbital wave theory.

###  IV.1 Parton mean-field theory

For small $K_{q}^{\prime}$ and $K_{o}^{\prime}$, we can start from the
solution in Sec. III to investigate quantum phase transitions out of the MSL
state. We first notice that the quartic terms contained in $H_{I}$ involve
only the matter fermions $\theta^{\gamma}$. Thus, at least at weak coupling
they do not alter the gauge-flux configuration in the ground state, and we can
work in the sector with fixed $W_{p}=1$ for all plaquettes.

Figure 5: Schematic representation of the ferroquadrupolar-vortex state on the
honeycomb lattice. The arrows represent the direction of the vector defined by
$(\langle O_{2}\rangle,\langle
O_{3}\rangle)=2(\langle\tau^{x}\rangle,\langle\tau^{z}\rangle)\propto(\cos\phi_{n},-\sin\phi_{n})$,
where $\phi_{n}=(n-1)\pi/3$, with $n=1,\dots,6$ being the sublattice index.

Next, we treat the interactions within a parton mean-field approach Savary and
Balents (2017); Jin  _et al._ (2022). To search for multipolar orders, we use
the mean-field decoupling

$\displaystyle\tilde{O}^{\gamma}_{j}\tilde{O}^{\gamma}_{l}$
$\displaystyle\mapsto$
$\displaystyle\langle\tilde{O}^{\gamma}_{j}\rangle\tilde{O}^{\gamma}_{l}+\langle\tilde{O}^{\gamma}_{l}\rangle\tilde{O}^{\gamma}_{j}-\langle\tilde{O}^{\gamma}_{j}\rangle\langle\tilde{O}^{\gamma}_{l}\rangle,$
(20) $\displaystyle T^{xyz}_{j}T^{xyz}_{l}$ $\displaystyle\mapsto$
$\displaystyle\langle T^{xyz}_{j}\rangle T^{xyz}_{l}+T^{xyz}_{j}\langle
T^{xyz}_{l}\rangle$ (21) $\displaystyle-\langle T^{xyz}_{j}\rangle\langle
T^{xyz}_{l}\rangle.$

The quadrupolar interaction $K_{q}^{\prime}$ is equivalent to an orbital
compass model on the honeycomb lattice Wu (2008); Khaliullin  _et al._ (2021)
and favors a six-sublattice FQV state that breaks rotational and translational
symmetries (see Fig. 5). We parametrize the quadrupolar order parameter as

$\langle\tilde{O}^{\gamma}_{j}\rangle=\rho\sum^{2}_{\mu=1}f^{\gamma}_{\mu}(\varphi)\cos(\mathbf{Q}_{\mu}\cdot\mathbf{R}_{j}).$
(22)

Here $\rho$ and $\varphi$ are variational parameters, the functions
$f^{\gamma}_{\mu}(\varphi)$ are given by
$f^{\gamma}_{1}(\varphi)=\sin(2\pi\gamma/3-\varphi)$ and
$f^{\gamma}_{2}(\varphi)=\cos(2\pi\gamma/3-\varphi)$, and the wave vectors
$\mathbf{Q}_{1}=\frac{\pi}{\sqrt{3}}(\sqrt{3},1)$ and
$\mathbf{Q}_{2}=\frac{\pi}{\sqrt{3}}(-\sqrt{3},1)$ connect high-symmetry
points of the reciprocal lattice Fernandes  _et al._ (2019). On the other
hand, for $K_{o}^{\prime}>0$ the octupolar order parameter $\langle
T^{xyz}_{j}\rangle$ is constrained to having the familiar form for
antiferromagnetic order on a bipartite lattice, i.e.,

$\langle T^{xyz}_{j}\rangle=\begin{dcases}+\chi\;\text{if $j\in\text{A}$},\\\
-\chi\;\text{if $j\in\text{B}$}.\end{dcases}$ (23)

A nonzero value of $\chi$ implies a phase with AO order.

Figure 6: Ground-state phase diagram obtained from the parton mean-field
theory setting $K_{o}=0.2K_{q}$. Here AO and FQV refer to the antiferro-
octupolar and ferroquadrupolar-vortex phases, respectively. In addition, the
system exhibits a phase with coexisting AO and FQV orders.

The precise form of the mean-field Hamiltonian can be found in Appendix B.
From this Hamiltonian, we obtain the free energy

$\displaystyle\mathcal{F}_{\text{MF}}$ $\displaystyle=$
$\displaystyle-T\sum^{18}_{\ell=1}\int_{\frac{1}{6}(\text{BZ})}\frac{d^{2}\mathbf{k}}{\mathcal{A}_{\text{BZ}}}\ln\big{[}1+e^{-\beta
E_{\ell}(\mathbf{k})}\big{]}$ (24)
$\displaystyle+K_{q}^{\prime}\big{[}\mathcal{A}_{0}-\mathcal{B}_{0}\sin(2\varphi)\big{]}\rho^{2}+3K_{o}^{\prime}\chi^{2},$

where $\beta=1/T$ is the inverse temperature and

$\displaystyle\mathcal{A}_{0}$
$\displaystyle\equiv\bigg{|}2\cos\bigg{(}\frac{\pi}{\sqrt{3}}\bigg{)}+\cos\bigg{(}\frac{2\pi}{\sqrt{3}}\bigg{)}\bigg{|}\approx
1.365,$ (25) $\displaystyle\mathcal{B}_{0}$
$\displaystyle\equiv\sqrt{3}\sin\bigg{(}\frac{\pi}{2\sqrt{3}}\bigg{)}\sin\bigg{(}\frac{\sqrt{3}\pi}{2}\bigg{)}\approx
0.557.$ (26)

In addition, $\mathcal{A}_{\text{BZ}}=8\pi^{2}/(3\sqrt{3})$ refers to the area
of the BZ and $E_{\ell}(\mathbf{k})$ are the dispersion relations obtained by
numerical diagonalization of the mean-field Hamiltonian.

Figure 7: Energy dispersion of the mean-field Hamiltonian for
$K_{o}=0.2K_{q}$, $K_{q}^{\prime}=1.3K_{q}$, and three different values of
$K_{o}^{\prime}$ in the broken-symmetry phases: (a) $K_{o}^{\prime}=0.65K_{q}$
in the AO phase, (b) $K_{o}^{\prime}=0.10K_{q}$ in the FQV phase, and (c)
$K_{o}^{\prime}=0.40K_{q}$ in the coexistence region. Here the Brillouin zone
is defined with respect to the unit cell shown in Fig. 1.

To determine the mean-field parameters $\rho$, $\varphi$, and $\chi$, we
numerically minimize the free energy in the zero-temperature limit by means of
a random-search algorithm which finds the lowest value of
$\mathcal{F}_{\text{MF}}$ from a set of randomly chosen initial conditions.
The ground-state phase diagram as a function of the integrability-breaking
perturbations for fixed $K_{o}=0.2K_{q}$ is shown in Fig. 6. First, we note
that the MSL phase is stable as long as $K_{q}$ and $K_{o}$ are both nonzero
and there are no zero-energy flat bands in the spectrum of the unperturbed
model. As we increase $K^{\prime}_{o}$, the MSL undergoes a quantum phase
transition to a state with AO order, $\chi\neq 0$, which spontaneously breaks
time-reversal symmetry. On the other hand, if we start from the solvable model
and increase $K^{\prime}_{q}$, the system develops FQV order with $\rho\neq
0$, breaking lattice rotation and translation symmetries. In this regime, we
find that the parameter $\varphi$ evaluates to either $\varphi=\pi/4$ or
$\varphi=5\pi/4$. This degeneracy is associated with the invariance of the
Hamiltonian under inversion of all quadrupole moments, which follows from
$f^{\gamma}_{\mu}(\varphi+\pi)=-f^{\gamma}_{\mu}(\varphi)$. Finally, when
$K_{q}^{\prime}$ and $K_{o}^{\prime}$ are comparable, the phase diagram
exhibits a coexistence region between the FQV and AO phases. Within our mean-
field approach, all phase transitions are of first order.

We now analyze the dispersion relation of the Majorana fermions in the ordered
phases. Representative results for the AO, FQV, and coexistence phases are
shown in Fig. 7. In the AO phase, the mean-field decoupling with $\chi\neq 0$
introduces a staggered on-site term that couples $\theta^{x}$ and $\theta^{z}$
without affecting $\theta^{y}$. As a result, the spectrum remains gapless with
a Dirac node at the $\mathbf{K}$ point [see Fig. 7(a)]. By contrast, in the
FQV phase the mean-field Hamiltonian with $\rho\neq 0$ couples all three
flavors of Majorana fermions. In this case, we obtain a fully gapped spectrum
[see Fig. 7(b)]. In the coexistence region, we find that the gap closes as one
of the bands crosses zero energy, forming a Majorana Fermi surface [see Fig.
7(c)]. Note that the spectrum is not particle-hole symmetric because the on-
site terms break the chiral symmetry in the mean-field Hamiltonian.

###  IV.2 Orbital wave theory

In the regime $K_{q}^{\prime},K_{o}^{\prime}\gg K_{q},K_{o}$, the mean-field
approximation used in the previous section is no longer reliable because the
strong integrability-breaking interactions are able to create vison
excitations. However, we can go the other way and consider the MSL Hamiltonian
$H_{s}$ to be a small perturbation to $H_{I}$. It is instructive to rewrite
the Hamiltonian in terms of the pseudospin and pseudo-orbital operators using
Eq. (8):

$\displaystyle H_{I}$ $\displaystyle=$ $\displaystyle
4K_{q}^{\prime}\sum_{\gamma}\sum_{\langle
jl\rangle_{\gamma}}(\mathbf{u}^{\gamma}\cdot\bm{\tau}_{j})(\mathbf{u}^{\gamma}\cdot\bm{\tau}_{l})+4K_{o}^{\prime}\sum_{\langle
jl\rangle}\tau^{y}_{j}\tau^{y}_{l},$ $\displaystyle H_{s}$ $\displaystyle=$
$\displaystyle 16K_{q}\sum_{\gamma}\sum_{\langle
jl\rangle_{\gamma}}(s^{\gamma}_{j}\tau^{y}_{j})(s^{\gamma}_{l}\tau^{y}_{l})$
(27) $\displaystyle+16K_{o}\sum_{\gamma}\sum_{\langle
jl\rangle_{\gamma}}(s^{\gamma}_{j}\mathbf{v}^{\gamma}\cdot\bm{\tau}_{j})(s^{\gamma}_{l}\mathbf{v}^{\gamma}\cdot\bm{\tau}_{l}),$

where
$\mathbf{u}^{\gamma}={\mathbf{v}}^{\gamma}\times\hat{\mathbf{y}}=\sin(2\pi\gamma/3)\hat{\mathbf{x}}+\cos(2\pi\gamma/3)\hat{\mathbf{z}}$.
Clearly, $H_{I}$ describes interactions within the pseudo-orbital sector. For
$K_{q}^{\prime},K_{o}^{\prime}\gg K_{q},K_{o}$, we expect $\bm{\tau}_{j}$ to
develop long-range order first, and the resulting ordered state determines the
effective interactions in the pseudospin sector at a lower energy scale.

Let us consider the classical ground states of $H_{I}$, with pseudo-orbitals
treated as vectors with length $\tau$. The FQV phase corresponds to a vortex
state with $\bm{\tau}$ vectors in the $xz$ plane, as illustrated in Fig. 5. In
the AO phase, the $\bm{\tau}$ vectors point along the $y$ axis. If we consider
a state that interpolates between FQV and AO configurations, with orbital
vectors forming an angle $\theta$ with the $\pm\hat{\mathbf{y}}$ direction
(alternating between A and B sublattices), the corresponding energy per site
is

$E_{\rm
cl}(\theta)=-3\tau^{2}(K_{q}^{\prime}\sin^{2}\theta+2K_{o}^{\prime}\cos^{2}\theta).$
(28)

Minimizing this energy, we obtain a direct transition between AO order
($\theta=0$) and FQV order ($\theta=\pi/2$) at
$K_{o}^{\prime}=K^{\prime}_{q}/2$. In contrast to the mean-field theory in
Sec. IV.1, this classical analysis in the strong-coupling limit does not
support the presence of a coexistence region in the phase diagram.

We can go beyond the classical analysis and compute quantum corrections using
linear orbital wave theory Murthy  _et al._ (1997); Wu (2008); Zhao and Liu
(2008). The first step is to rotate the pseudo-orbital operators to the local
polarization axis in either FQV or AO states. We write
$\bm{\tau}_{j}=R^{T}_{j}\tilde{\bm{\tau}}_{j}$, where $R_{j}$ is a $3\times 3$
orthogonal matrix that implements the appropriate position-dependent rotation.
In the rotated frame, we employ the Holstein-Primakoff transformation

$\displaystyle\tilde{\tau}^{z}_{j}$ $\displaystyle=$
$\displaystyle\tau-a^{\dagger}_{j}a^{\phantom{\dagger}}_{j},$ (29)
$\displaystyle\tilde{\tau}^{-}_{j}$ $\displaystyle=$ $\displaystyle
a^{\dagger}_{j}(2\tau-a^{\dagger}_{j}a^{\phantom{\dagger}}_{j})^{1/2},$ (30)

where $a^{\dagger}_{j}$ and $a^{\phantom{\dagger}}_{j}$ are bosonic creation
and annihilation operators. The effective Hamiltonian to order $\tau$ is
quadratic in the bosonic operators and can be diagonalized by a Bogoliubov
transformation. We then obtain

$H_{I}\approx\sum_{\nu,\mathbf{k}}\omega_{\nu}(\mathbf{k})\alpha_{\nu}^{\dagger}(\mathbf{k})\alpha_{\nu}^{\phantom{\dagger}}(\mathbf{k})+E_{0},$
(31)

where $\omega_{\nu}(\mathbf{k})$ is the dispersion relation of the quantized
pseudo-orbital waves in the FQV or AO states, with $\nu$ being the band index
and $\alpha_{\nu}(\mathbf{k})$ being the corresponding boson annihilation
operator, and

$E_{0}=-3N\tau(\tau+1)\mathcal{E}+\frac{1}{2}\sum_{\nu,\mathbf{k}}\omega_{\nu}(\mathbf{k})$
(32)

is the ground state energy including the quantum correction, with
$\mathcal{E}=K_{q}^{\prime}$ in the FQV state and
$\mathcal{E}=2K_{o}^{\prime}$ in the AO state.

Perturbative stability of the ordered states against quantum fluctuations
requires $\omega_{\nu}(\mathbf{k})$ to be real and positive for all bands.
Using this criterion, we find that the FQV state remains locally stable for
$K_{o}^{\prime}<0.746K_{q}^{\prime}$, beyond the point
$K_{o}^{\prime}=K_{q}^{\prime}/2$ where its classical energy becomes higher
than that of the AO state. Consistently, the quantum correction to the ground-
state energy of the AO state is real for $K_{o}^{\prime}>K_{q}^{\prime}/2$
(see Fig. 8). In the AO phase, the pseudo-orbital wave spectrum is gapped, as
expected since the ground state breaks only discrete symmetries. In the FQV
phase, the spectrum contains a zero-energy flat band. As emphasized in Ref. Wu
(2008), the classical FQV state has an SO(2) degeneracy which is lifted by
quantum fluctuations, and the quantum corrections select states that maximize
the number of zero-energy modes. The gap in the low-energy band in the FQV
phase appears only beyond the linear orbital wave theory Khaliullin  _et al._
(2021). Since the excitation spectrum indicates a perturbative stability of
the ordered states in the corresponding regimes, with an energy level crossing
around $K_{o}^{\prime}=K_{q}^{\prime}/2$, we conclude that the orbital wave
theory points to a first-order transition between FQV and AO phases.

Figure 8: Quantum correction to the ground-state energy of FQV and AO states
[see Eqs. (28) and (32)]. The dispersion relation of pseudo-orbital waves in
the FQV state becomes complex for $K_{o}^{\prime}>0.746K_{q}^{\prime}$,
signaling an instability. Similarly, the AO state becomes unstable for
$K_{o}^{\prime}<0.5K_{q}^{\prime}$.

We can now derive an effective Hamiltonian for the pseudospins using a mean-
field decoupling of $H_{s}$ in which the pseudo-orbital degrees of freedom are
frozen in the ground state of $H_{I}$. This is a standard approach to analyze
the low-energy excitations of spin-orbital models van den Brink  _et al._
(1998). We obtain from Eq. (27)

$H^{\text{eff}}_{s}=\sum_{\gamma=x,y,z}\sum_{\langle
jl\rangle_{\gamma}}\mathcal{K}^{\gamma}_{jl}s^{\gamma}_{j}s^{\gamma}_{l},$
(33)

where

$\mathcal{K}^{\gamma}_{jl}=-16K_{o}\langle\tau^{y}_{j}\tau^{y}_{l}\rangle-16K_{q}\langle(\mathbf{v}^{\gamma}\cdot\bm{\tau}_{j})(\mathbf{v}^{\gamma}\cdot\bm{\tau}_{l})\rangle.$
(34)

The effective Hamiltonian in Eq. (33) has the precise form of a Kitaev model
for pseudospin-$\frac{1}{2}$ moments. As a consequence, the fate of the system
in the regime $K_{q}^{\prime},K_{o}^{\prime}\gg K_{q},K_{o}$ is to exhibit a
_hidden_ Kitaev spin liquid on top of classical AO or FQV phases. This feature
bears some resemblance to the behavior of the spin-$\frac{3}{2}$ Kitaev model
studied in Ref. Jin  _et al._ (2022), in which a spin-$\frac{1}{2}$ Kitaev
model emerges in the presence of a single-ion anisotropy that induces a
uniform spin quadrupole moment. In our case, the Kitaev spin liquid arises
either from octupolar order for $K^{\prime}_{o}\gg K_{q}^{\prime}$ or from
translation-symmetry-breaking quadrupolar order for $K^{\prime}_{q}\gg
K_{o}^{\prime}$.

The pseudospin excitation spectrum depends on the effective Kitaev couplings
in Eq. (34). We determine these couplings by computing the pseudo-orbital
correlations within the linear orbital wave theory. In the AO phase, which
preserves the $\mathbb{Z}_{3}$ symmetry, we obtain homogeneous
antiferromagnetic coupling $\mathcal{K}^{\gamma}_{jl}=\mathcal{K}>0$. As a
consequence, in the AO phase the Kitaev spin liquid harbors gapless Majorana
fermion excitations with linear dispersion about the $\mathbf{K}$ point. This
conclusion agrees with the result from the parton mean-field approach in Sec.
IV.1. Remarkably, the Kitaev spin liquid remains gapless despite the
spontaneous breakdown of time-reversal symmetry. In contrast to the effect of
a magnetic field in the original Kitaev model Kitaev (2006), here the AO order
parameter preserves the symmetry combining time reversal with a $C_{2}$
rotation that exchanges A and B sublattices. This symmetry rules out the
three-spin interaction that would drive a nonzero spin chirality in the Kitaev
spin liquid and generate a mass term for the Majorana fermions Kitaev (2006).

In the FQV phase, the effective Kitaev couplings $\mathcal{K}_{jl}^{\gamma}$
become bond dependent. While the FQV order breaks the $\mathbb{Z}_{3}$
symmetry, it preserves a $\mathbb{Z}_{6}$ symmetry corresponding to a rotation
in real and internal space around the center of the hexagon that defines the
enlarged unit cell in Fig. 5. This symmetry implies that there are only two
types of Kitaev couplings: $\mathcal{K}^{\gamma}_{jl}=\mathcal{K}_{\rm in}$
for bonds inside the unit cell, represented by solid lines in Fig. 5, and
$\mathcal{K}^{\gamma}_{jl}=\mathcal{K}_{\rm out}$ for bonds connecting
neighboring unit cells, represented by dashed lines. Using the linear orbital
wave theory, we find $\mathcal{K}_{\rm in}>0$ and $\mathcal{K}_{\rm out}>0$
with $\mathcal{K}_{\rm out}\ll\mathcal{K}_{\rm in}$ regardless of the ratio
between $K_{o}$ and $K_{q}$. In this regime, the Majorana spectrum is gapped
with weakly dispersive bands. This gapped FQV phase is consistent with the
result from the parton mean-field approach. We note that the effect of FQV
order on the Majorana fermions in the hidden Kitaev spin liquid is analogous
to a Kekulé distortion in graphene, which also triples the size of the unit
cell and opens a gap in the spectrum Hou  _et al._ (2007); A. C. Qu and P.
Nigge and S. Link and G. Levy and M. Michiardi and P. L. Spandar and T. Matthé
and M. Schneider and S. Zhdanovich and U. Starke and C. Gutiérrez and A.
Damascelli (2022).

##  V Conclusions

In summary, we investigated the properties of a model for
$j_{\mathrm{eff}}=\frac{3}{2}$ local moments on a honeycomb lattice with bond-
directional quadrupolar and octupolar interactions. For special values of the
parameters, we obtained an exact solution by representing the
$j_{\mathrm{eff}}=\frac{3}{2}$ multipolar operators in terms of two sets of
Majorana fermions. This representation allowed us to map the model onto
noninteracting fermions hopping in the background of a static $\mathbb{Z}_{2}$
gauge field. This model provides a prime example of an interacting system with
a nondegenerate ground state, gapless fermionic excitations, and QSL
correlations due to the fractionalization of higher-order multipoles.

We demonstrated the emergence of nontrivial topological phases by explicitly
breaking rotational and time-reversal symmetries in the MSL. Applying both
strain and magnetic fields, we found Abelian and non-Abelian phases with Chern
numbers of $0$, $\pm 1$, and $\pm 2$. The Abelian topological phase with
$\mathcal{N}_{\mathrm{Ch}}=\pm 2$ appears for arbitrarily weak perturbations.
While Abelian and non-Abelian phases were reported in recent studies of the
antiferromagnetic Kitaev model in a magnetic field Jiang  _et al._ (2020);
Zhang  _et al._ (2022), our model highlights the role of strain as an
additional control parameter for driving topological phase transitions between
chiral spin-orbital liquid phases.

The stability of the MSL was probed by considering integrability-breaking
quadrupolar and octupolar interactions. In the weak-coupling regime, the
application of parton mean-field theory showed that the MSL remains stable
over a wide range of interaction strengths. This stability can be understood
from a renormalization-group point of view by noting that in the fermionic
representation the interactions are clearly irrelevant due to the vanishing
density of states of the Dirac nodes in the MSL. When a phase transition out
of the MSL finally takes place, the system develops either ferroquadrupolar-
vortex or antiferro-octupolar order. The parton mean-field theory leads to a
coexistence region at intermediate couplings, but the orbital wave theory at
strong coupling indicates a direct first-order transition between the two
ordered phases. Despite the spontaneous symmetry breaking, the low-energy
spectrum in the strong-coupling regime can be interpreted in terms of a hidden
Kitaev spin liquid with effective couplings between pseudospin degrees of
freedom. Of course, an important open question is whether the regime of
dominant quadrupolar and octupolar interactions considered in this work can be
reached in layered materials with $4d$ and $5d$ transition-metal ions Takayama
_et al._ (2021).

###### Acknowledgements.

We would like to thank C. S. de Farias, E. Miranda, and W. Natori for
interesting discussions with. H.F. acknowledges funding from the Brazilian
agency CNPq under Grant No. 311428/2021-5. R.G.P. acknowledges funding from
CNPq under Grant No. 309569/2022-2. This work was supported by a grant from
the Simons Foundation (Grant No. 1023171, R.G.P.). Research at IIP-UFRN is
supported by the Brazilian ministries Ministério da Educação (MEC) and
Ministério da Ciência, Tecnologia e Inovação (MCTI).

##  APPENDIX A COMPARISON BETWEEN THE ENERGIES OF THE ZERO- AND $\pi$-GAUGE
FLUX GROUND STATES

In the main text, we stated that the ground state (GS) of the MSL lies in the
zero-gauge flux sector $W_{p}=+1$, which is obtained by setting $u_{\langle
jl\rangle_{\gamma}}=1$ for sites $j$ in sublattice A and $l$ in sublattice B.
To confirm that result, we will compare the GS energy of the MSL for both
$W_{p}=+1$ and $W_{p}=-1$, in which the latter defines the so-called
$\pi$-gauge flux sector for this model.

Figure 9: (a) Honeycomb lattice for a translationally invariant $\pi$-gauge
flux configuration $W_{p}=-1$. The dashed blue and filled red bonds refer to
the bond variables $u_{\langle jl\rangle_{\gamma}}=+1$ and $u_{\langle
jl\rangle_{\gamma}}=-1$, respectively. In terms of this configuration, the
physical unit cell contains four sublattices, which are indicated by numbers
and also by different regular polygons. (b) Behavior of the MSL GS energies
$E^{+}_{\mathrm{GS}}\equiv E_{\mathrm{GS}}(\\{W_{p}=+1\\})$ and
$E^{-}_{\mathrm{GS}}\equiv E_{\mathrm{GS}}(\\{W_{p}=-1\\})$ for the zero- and
$\pi$-gauge flux configurations, respectively. Note that $E^{+}_{\mathrm{GS}}$
is always lower than $E^{-}_{\mathrm{GS}}$, which is consistent with Lieb’s
theorem Lieb (1994).

As exemplified in Fig. 9(a), a translationally invariant honeycomb lattice
with $W_{p}=-1$ is obtained with a unit cell containing four sublattice sites.
Since the gauge variables for this unit cell evaluate to

$\displaystyle u_{\langle 3,2\rangle_{x}}=+1,\hskip 28.45274ptu_{\langle
3,2\rangle_{y}}=+1,\hskip 28.45274ptu_{\langle 1,2\rangle_{z}}=+1,$ (35)
$\displaystyle u_{\langle 1,4\rangle_{x}}=+1,\hskip 28.45274ptu_{\langle
1,4\rangle_{y}}=-1,\hskip 28.45274ptu_{\langle 3,4\rangle_{z}}=+1,$ (36)

we find that the MSL Hamiltonian in Eq. (14) turns out to be

$\displaystyle H_{s}=$
$\displaystyle\sum_{\mathbf{k}\in\overline{\mathrm{BZ}}}\bar{\bm{\Theta}}^{\dagger}(\mathbf{k})\begin{pmatrix}0&i\Gamma^{z}(\mathbf{k})&0&i[\Gamma^{x}(\mathbf{k})-\Gamma^{y}(\mathbf{k})]\\\
-i\Gamma^{z\dagger}(\mathbf{k})&0&-i[\Gamma^{x\dagger}(\mathbf{k})+\Gamma^{y\dagger}(\mathbf{k})]&0\\\
0&i[\Gamma^{x}(\mathbf{k})+\Gamma^{y}(\mathbf{k})]&0&i\Gamma^{z}(\mathbf{k})\\\
-i[\Gamma^{x\dagger}(\mathbf{k})-\Gamma^{y\dagger}(\mathbf{k})]&0&-i\Gamma^{z\dagger}(\mathbf{k})&0\end{pmatrix}\bar{\bm{\Theta}}(\mathbf{k}),$
(37)

where
$\bar{\bm{\Theta}}(\mathbf{k})\equiv\begin{pmatrix}\bm{\theta}_{1}(\mathbf{k}),&\bm{\theta}_{2}(\mathbf{k}),&\bm{\theta}_{3}(\mathbf{k}),&\bm{\theta}_{4}(\mathbf{k})\end{pmatrix}^{T}$;
$\overline{\mathrm{BZ}}$ refers to the Brillouin zone associated with the
four-sublattice unit cell; and the functions $\Gamma^{\gamma}(\mathbf{k})$,
with $\gamma\in\\{x,y,z\\}$, represent the following momentum-dependent
matrices:

$\displaystyle\Gamma^{x}(\mathbf{k})=\Bigg{[}K_{q}\begin{pmatrix}0&0&0\\\
0&1&0\\\ 0&0&0\end{pmatrix}+\frac{K_{o}}{4}\begin{pmatrix}1&0&\sqrt{3}\\\
0&0&0\\\
\sqrt{3}&0&3\end{pmatrix}\Bigg{]}e^{-i\mathbf{k}\cdot\bm{\delta}_{1}},$ (38)
$\displaystyle\Gamma^{y}(\mathbf{k})=\Bigg{[}K_{q}\begin{pmatrix}0&0&0\\\
0&1&0\\\ 0&0&0\end{pmatrix}+\frac{K_{o}}{4}\begin{pmatrix}1&0&-\sqrt{3}\\\
0&0&0\\\
-\sqrt{3}&0&3\end{pmatrix}\Bigg{]}e^{-i\mathbf{k}\cdot\bm{\delta}_{2}},$ (39)
$\displaystyle\Gamma^{z}(\mathbf{k})=\Bigg{[}K_{q}\begin{pmatrix}0&0&0\\\
0&1&0\\\ 0&0&0\end{pmatrix}+K_{o}\begin{pmatrix}1&0&0\\\ 0&0&0\\\
0&0&0\end{pmatrix}\Bigg{]}e^{-i\mathbf{k}\cdot\bm{\delta}_{3}}.$ (40)

Here, $\bm{\delta}_{1}$, $\bm{\delta}_{2}$, and $\bm{\delta}_{3}$ are nearest-
neighbor vectors pointing from an even- to an odd-sublattice site in the
direction of the $x$, $y$, and $z$ bonds, respectively.

The diagonalization of $H_{s}$ in Eq. (37) yields $3\bar{N}_{s}=12$ energy
dispersions $\bar{E}_{s,n}(\mathbf{k})$, but only four of them can be computed
in closed form. Since the latter also have a complicated dependence, they will
not be displayed here. Regardless, the GS energy $E^{-}_{\mathrm{GS}}\equiv
E_{\mathrm{GS}}(\\{W_{p}=-1\\})$ associated with this Hamiltonian can be
computed numerically. In fact, its formula is given by

$\frac{E^{-}_{\mathrm{GS}}}{N}=\frac{1}{\bar{N}_{s}}\sum\limits^{3\bar{N}_{s}}_{n=1}\int_{\overline{\mathrm{BZ}}}\frac{d^{2}\mathbf{k}}{\bar{\mathcal{A}}_{\overline{\mathrm{BZ}}}}\bar{E}_{n}(\mathbf{k})\varTheta[-\bar{E}_{s,n}(\mathbf{k})],$
(41)

where $N$ refers to the number of sites of the honeycomb lattice,
$\bar{\mathcal{A}}_{\overline{\mathrm{BZ}}}=4\pi^{2}/(3\sqrt{3}a^{2})$ is the
area of the $\overline{\mathrm{BZ}}$, and $\varTheta(x)$ denotes the Heaviside
step function. This formula should be contrasted with the equation

$\frac{E^{+}_{\mathrm{GS}}}{N}=\frac{1}{N_{s}}\sum\limits^{3N_{s}}_{n=1}\int_{\mathrm{BZ}}\frac{d^{2}\mathbf{k}}{\mathcal{A}_{\mathrm{BZ}}}E_{n}(\mathbf{k})\varTheta[-E_{s,n}(\mathbf{k})],$
(42)

which yields the GS energy $E^{+}_{\mathrm{GS}}\equiv
E_{\mathrm{GS}}(\\{W_{p}=+1\\})$ of the MSL in the zero-gauge flux
configuration. In this case, the unit cell has $N_{s}=2$ sublattices, and
$\mathcal{A}_{\mathrm{BZ}}=8\pi^{2}/(3\sqrt{3}a^{2})$ gives the area of the
corresponding BZ. In addition, $E_{s,n}(\mathbf{k})$ refers to $3N_{s}=6$
energy dispersions whose expressions can be found in the main text [see Eq.
(16) and the paragraph above it].

In Fig. 9(b), we show the behavior of both $E^{-}_{\mathrm{GS}}$ and
$E^{+}_{\mathrm{GS}}$ as a function of the ratio $|K_{o}/K_{q}|$. Clearly, the
zero-gauge flux configuration $W_{p}=+1$ has the lowest GS energy. Moreover,
we also investigate the behavior of the GS energy gap per lattice site
$(E^{-}_{\mathrm{GS}}-E^{+}_{\mathrm{GS}})/N$. In fact, our numerical results
reveal that $(E^{-}_{\mathrm{GS}}-E^{+}_{\mathrm{GS}})/N$ depends linearly on
$|K_{o}/K_{q}|$. Consequently, we obtain

$\frac{E^{-}_{\mathrm{GS}}-E^{+}_{\mathrm{GS}}}{N}=b_{0}(|K_{o}|+2|K_{q}|),$
(43)

where $b_{0}\approx 0.012867$. Remarkably, this means that the quadrupolar and
octupolar degrees of freedom of the MSL contribute differently to
$E^{-}_{\mathrm{GS}}-E^{+}_{\mathrm{GS}}$. In other words, the impact of the
former is twice that of the latter.

Note that, according to Lieb’s theorem Lieb (1994), the zero-gauge flux
configuration yields the lowest-energy GS for mirror-symmetric models with
nearest-neighbor interactions defined on a honeycomb lattice. In light of this
result, consider the MSL Hamiltonian $H_{s}$ in Eq. (14). The first term in
$H_{s}$ involving the $\theta^{y}$ fermions is clearly invariant under a
reflection exchanging the $x$ and $y$ bonds since $\theta^{y}$ transforms as
the matter Majorana fermions in the Kitaev model for isotropic exchange
interactions. However, the mirror symmetry of the term that involves
$\theta^{x}$ and $\theta^{z}$ requires the additional transformation
$\theta^{z}\to-\theta^{z}$, which can be traced back to the symmetry
properties of the quadrupole components $(\tau^{x},\;\tau^{z})$. We can view
$\theta^{z}\to-\theta^{z}$ as a gauge transformation which does not affect the
gauge-invariant fluxes nor the spectrum of $H_{s}$. For this reason, we
expected the result of Lieb’s theorem to apply also to our model, which was
confirmed by our numerical results.

##  APPENDIX B MEAN-FIELD HAMILTONIAN FOR MAJORANA FERMIONS IN THE ZERO-GAUGE
FLUX SECTOR

After we fix the uniform gauge configuration for $u_{\langle jl\rangle}$,
which is equivalent to setting $W_{p}=+1$ for all plaquettes, the exactly
solvable Hamiltonian can be written as

$H_{s}+\delta
H_{s}=\sum_{\mathbf{k}}\bm{\Theta}^{\dagger}(\mathbf{k})\bm{\mathcal{H}}_{s}(\mathbf{k})\bm{\Theta}(\mathbf{k}),$
(44)

where
$\bm{\Theta}(\mathbf{k})\equiv(\theta^{x}_{A}(\mathbf{k}),\theta^{x}_{B}(\mathbf{k}),\theta^{y}_{A}(\mathbf{k}),\theta^{y}_{B}(\mathbf{k}),\theta^{z}_{A}(\mathbf{k}),\theta^{z}_{B}(\mathbf{k}))^{T}$
and we define the matrix

$\displaystyle\bm{\mathcal{H}}_{s}(\mathbf{k})=\begin{pmatrix}0&iK_{o}g_{x}(\mathbf{k})&i\varepsilon_{3}&0&-ih&iK_{o}g_{xz}(\mathbf{k})\\\
-K_{o}g^{*}_{x}(\mathbf{k})&0&0&i\varepsilon_{3}&-iK_{o}g^{*}_{xz}(\mathbf{k})&-ih\\\
-i\varepsilon_{3}&0&0&iK_{q}g_{y}(\mathbf{k})&i\varepsilon_{2}&0\\\
0&-i\varepsilon_{3}&-iK_{q}g^{*}_{y}(\mathbf{k})&0&0&i\varepsilon_{2}\\\
ih&iK_{o}g_{xz}(\mathbf{k})&-i\varepsilon_{2}&0&0&iK_{o}g_{z}(\mathbf{k})\\\
-iK_{o}g^{*}_{xz}(\mathbf{k})&ih&0&-i\varepsilon_{2}&-iK_{o}g_{z}(\mathbf{k})&0\end{pmatrix},$
(45)

with the functions given below Eq. (16). Here we have included the
perturbations that determine the phase diagram in Fig. 3. In the parton mean-
field theory in Sec. IV.1, we set $\varepsilon_{2}=\varepsilon_{3}=h=0$.

According to Eqs. (22) and (23), the mean-field decoupling of the Hamiltonian
$H=H_{s}+H_{I}$ leads to

$\displaystyle H_{\text{MF}}=$
$\displaystyle\sum_{\mathbf{k}\in\frac{1}{6}(\text{BZ})}\begin{pmatrix}\bm{\Theta}(\mathbf{k})\\\
\bm{\Theta}(\mathbf{k}+\mathbf{Q}_{1})\\\
\bm{\Theta}(\mathbf{k}+\mathbf{Q}_{2})\end{pmatrix}^{\dagger}\begin{pmatrix}\bm{\mathcal{H}}_{s}(\mathbf{k})+\bm{\mathcal{V}}&\bm{\mathcal{V}}_{\mathbf{Q}_{1}}&\bm{\mathcal{V}}_{\mathbf{Q}_{2}}\\\
\bm{\mathcal{V}}^{\dagger}_{\mathbf{Q}_{1}}&\bm{\mathcal{H}}_{s}(\mathbf{k}+\mathbf{Q}_{1})+\bm{\mathcal{V}}&\mathbb{0}\\\
\bm{\mathcal{V}}^{\dagger}_{\mathbf{Q}_{2}}&\mathbb{0}&\bm{\mathcal{H}}_{s}(\mathbf{k}+\mathbf{Q}_{2})+\bm{\mathcal{V}}\end{pmatrix}\begin{pmatrix}\bm{\Theta}(\mathbf{k})\\\
\bm{\Theta}(\mathbf{k}+\mathbf{Q}_{1})\\\
\bm{\Theta}(\mathbf{k}+\mathbf{Q}_{2})\end{pmatrix}$
$\displaystyle-\frac{K_{q}^{\prime}}{2}N\rho^{2}\sum^{2}_{\mu=1}\sum_{\gamma=1}^{3}[f_{\mu}^{\gamma}(\varphi)]^{2}\cos(\mathbf{Q}_{\mu}\cdot\bm{\delta}_{\gamma})+\frac{3K_{o}^{\prime}}{2}N\chi^{2},$
(46)

where we have set

$\displaystyle\bm{\mathcal{V}}_{\mathbf{Q}_{\mu}}$
$\displaystyle=3K_{q}^{\prime}\begin{pmatrix}0&0&i\rho^{xy*}_{\mu,A}&0&0&0\\\
0&0&0&i\rho^{xy}_{\mu,B}&0&0\\\
-i\rho^{xy*}_{\mu,A}&0&0&0&-i\rho^{yz*}_{\mu,A}&0\\\
0&-i\rho^{xy}_{\mu,B}&0&0&0&-i\rho^{yz}_{\mu,B}\\\
0&0&i\rho^{yz*}_{\mu,A}&0&0&0\\\ 0&0&0&i\rho^{yz}_{\mu,B}&0&0\end{pmatrix},$
(47) $\displaystyle\bm{\mathcal{V}}$
$\displaystyle=6K_{o}^{\prime}\begin{pmatrix}0&0&0&0&-i\chi&0\\\
0&0&0&0&0&i\chi\\\ 0&0&0&0&0&0\\\ 0&0&0&0&0&0\\\ i\chi&0&0&0&0&0\\\
0&-i\chi&0&0&0&0\end{pmatrix}.$ (48)

The matrices $\mathcal{V}_{\mathbf{Q}_{\mu}}$ are defined in terms of the
functions

$\displaystyle\rho^{xy}_{\mu,A}$
$\displaystyle=\frac{\rho}{2}\left[\frac{f^{1}_{\mu}(\varphi)e^{i\mathbf{Q}_{\mu}\cdot\bm{\delta}_{1}}+f^{2}_{\mu}(\varphi)e^{i\mathbf{Q}_{\mu}\cdot\bm{\delta}_{2}}}{2}-f^{3}_{\mu}(\varphi)e^{i\mathbf{Q}_{\mu}\cdot\bm{\delta}_{3}}\right],$
(49) $\displaystyle\rho^{yz}_{\mu,A}$
$\displaystyle=\frac{\sqrt{3}\rho}{4}\left[-f^{1}_{\mu}(\varphi)e^{i\mathbf{Q}_{\mu}\cdot\bm{\delta}_{1}}+f^{2}_{\mu}(\varphi)e^{i\mathbf{Q}_{\mu}\cdot\bm{\delta}_{2}}\right],$
(50) $\displaystyle\rho^{xy}_{\mu,B}$
$\displaystyle=\frac{\rho}{2}\left[\frac{f^{1}_{\mu}(\varphi)e^{i\mathbf{Q}_{\mu}\cdot\mathbf{n}_{1}}+f^{2}_{\mu}(\varphi)e^{i\mathbf{Q}_{\mu}\cdot\mathbf{n}_{2}}}{2}-f^{3}_{\mu}(\varphi)\right],$
(51) $\displaystyle\rho^{yz}_{\mu,B}$
$\displaystyle=\frac{\sqrt{3}\rho}{4}\left[-f^{1}_{\mu}(\varphi)e^{i\mathbf{Q}_{\mu}\cdot\mathbf{n}_{1}}+f^{2}_{\mu}(\varphi)e^{i\mathbf{Q}_{\mu}\cdot\mathbf{n}_{2}}\right].$
(52)

## References

  * Balents (2010) L. Balents, Spin liquids in frustrated magnets, Nature 464, 199 (2010).
  * Broholm  _et al._ (2020) C. Broholm, R. J. Cava, S. A. Kivelson, D. G. Nocera, M. R. Norman, and T. Senthil, Quantum spin liquids, Science 367, eaay0668 (2020).
  * Savary and Balents (2017) L. Savary and L. Balents, Quantum spin liquids: a review, Rep. Prog. Phys. 80, 016502 (2017).
  * Takagi  _et al._ (2019) H. Takagi, T. Takayama, G. Jackeli, G. Khaliullin, and S. E. Nagler, Concept and realization of Kitaev quantum spin liquids, Nat. Rev. Phys. 1, 264 (2019).
  * Trebst and Hickey (2022) S. Trebst and C. Hickey, Kitaev materials, Phys. Rep. 950, 1 (2022).
  * Kitaev (2006) A. Kitaev, Anyons in an exactly solved model and beyond, Ann. Phys. 321, 2 (2006).
  * Jackeli and Khaliullin (2009) G. Jackeli and G. Khaliullin, Mott Insulators in the Strong Spin-Orbit Coupling Limit: From Heisenberg to a Quantum Compass and Kitaev Models, Phys. Rev. Lett. 102, 017205 (2009).
  * Singh and Gegenwart (2010) Y. Singh and P. Gegenwart, Antiferromagnetic Mott insulating state in single crystals of the honeycomb lattice material ${\text{Na}}_{2}{\text{IrO}}_{3}$, Phys. Rev. B 82, 064412 (2010).
  * Singh  _et al._ (2012) Y. Singh, S. Manni, J. Reuther, T. Berlijn, R. Thomale, W. Ku, S. Trebst, and P. Gegenwart, Relevance of the Heisenberg-Kitaev Model for the Honeycomb Lattice Iridates ${A}_{2}{\mathrm{IrO}}_{3}$, Phys. Rev. Lett. 108, 127203 (2012).
  * Cao  _et al._ (2013) G. Cao, T. F. Qi, L. Li, J. Terzic, V. S. Cao, S. J. Yuan, M. Tovar, G. Murthy, and R. K. Kaul, Evolution of magnetism in the single-crystal honeycomb iridates $({\text{Na}}_{1-x}$Li${}_{x}{)}_{2}\text{Ir}{\text{O}}_{3}$, Phys. Rev. B 88, 220414 (2013).
  * 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, A spin-orbital-entangled quantum liquid on a honeycomb lattice, Nature (London) 554, 341 (2018).
  * Plumb  _et al._ (2014) K. W. Plumb, J. P. Clancy, L. J. Sandilands, V. V. Shankar, Y. F. Hu, K. S. Burch, H.-Y. Kee, and Y.-J. Kim, $\alpha-{\mathrm{RuCl}}_{3}$: A spin-orbit assisted Mott insulator on a honeycomb lattice, Phys. Rev. B 90, 041112 (2014).
  * 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, Majorana quantization and half-integer thermal quantum Hall effect in a Kitaev spin liquid, Nature (London) 559, 227 (2018).
  * Tanaka  _et al._ (2022) O. Tanaka, Y. Mizukami, R. Harasawa, K. Hashimoto, K. Hwang, N. Kurita, H. Tanaka, S. Fujimoto, Y. Matsuda, E.-G. Moon, and T. Shibauchi, Thermodynamic evidence for a field-angle-dependent Majorana gap in a Kitaev spin liquid, Nat. Phys. 1, 1 (2022).
  * Witczak-Krempa  _et al._ (2014) W. Witczak-Krempa, G. Chen, Y. B. Kim, and L. Balents, Correlated Quantum Phenomena in the Strong Spin-Orbit Regime, Annu. Rev. Condens. Matter Phys. 5, 57 (2014).
  * Takayama  _et al._ (2021) T. Takayama, J. Chaloupka, A. Smerald, G. Khaliullin, and H. Takagi, Spin-Orbit-Entangled Electronic Phases in 4$d$ and 5$d$ Transition-Metal Compounds, J. Phys. Soc. Jpn. 90, 062001 (2021).
  * Stavropoulos  _et al._ (2019) P. P. Stavropoulos, D. Pereira, and H.-Y. Kee, Microscopic Mechanism for a Higher-Spin Kitaev Model, Phys. Rev. Lett. 123, 037203 (2019).
  * Xu  _et al._ (2020) C. Xu, J. Feng, M. Kawamura, Y. Yamaji, Y. Nahas, S. Prokhorenko, Y. Qi, H. Xiang, and L. Bellaiche, Possible Kitaev Quantum Spin Liquid State in 2D Materials with $S=3/2$, Phys. Rev. Lett. 124, 087205 (2020).
  * Xu  _et al._ (2018) C. Xu, J. Feng, H. Xiang, and L. Bellaiche, Interplay between Kitaev interaction and single ion anisotropy in ferromagnetic CrI3 and CrGeTe3 monolayers, npj Comput. Mater. 4, 57 (2018).
  * Stavropoulos  _et al._ (2021) P. P. Stavropoulos, X. Liu, and H.-Y. Kee, Magnetic anisotropy in spin-3/2 with heavy ligand in honeycomb Mott insulators: Application to ${\mathrm{CrI}}_{3}$, Phys. Rev. Res. 3, 013216 (2021).
  * Baskaran  _et al._ (2008) G. Baskaran, D. Sen, and R. Shankar, Spin-$S$ Kitaev model: Classical ground states, order from disorder, and exact correlation functions, Phys. Rev. B 78, 115116 (2008).
  * Koga  _et al._ (2018) A. Koga, H. Tomishige, and J. Nasu, Ground-state and Thermodynamic Properties of an $S=1$ Kitaev Model, J. Phys. Soc. Jpn. 87, 063703 (2018).
  * Hickey  _et al._ (2020) C. Hickey, C. Berke, P. P. Stavropoulos, H.-Y. Kee, and S. Trebst, Field-driven gapless spin liquid in the spin-1 Kitaev honeycomb model, Phys. Rev. Res. 2, 023361 (2020).
  * Dong and Sheng (2020) X.-Y. Dong and D. N. Sheng, Spin-1 Kitaev-Heisenberg model on a honeycomb lattice, Phys. Rev. B 102, 121102 (2020).
  * Zhu  _et al._ (2020) Z. Zhu, Z.-Y. Weng, and D. N. Sheng, Magnetic field induced spin liquids in $S=1$ Kitaev honeycomb model, Phys. Rev. Res. 2, 022047 (2020).
  * Khait  _et al._ (2021) I. Khait, P. P. Stavropoulos, H.-Y. Kee, and Y. B. Kim, Characterizing spin-one Kitaev quantum spin liquids, Phys. Rev. Res. 3, 013160 (2021).
  * Jin  _et al._ (2022) H.-K. Jin, W. M. H. Natori, F. Pollmann, and J. Knolle, Unveiling the $S=3/2$ Kitaev honeycomb spin liquids, Nat. Commun. 13, 3813 (2022).
  * Chen  _et al._ (2010) G. Chen, R. Pereira, and L. Balents, Exotic phases induced by strong spin-orbit coupling in ordered double perovskites, Phys. Rev. B 82, 174440 (2010).
  * Natori  _et al._ (2016) W. M. H. Natori, E. C. Andrade, E. Miranda, and R. G. Pereira, Chiral Spin-Orbital Liquids with Nodal Lines, Phys. Rev. Lett. 117, 017204 (2016).
  * Natori  _et al._ (2017) W. M. H. Natori, M. Daghofer, and R. G. Pereira, Dynamics of a $j=\frac{3}{2}$ quantum spin liquid, Phys. Rev. B 96, 125109 (2017).
  * Romhányi  _et al._ (2017) J. Romhányi, L. Balents, and G. Jackeli, Spin-Orbit Dimers and Noncollinear Phases in ${d}^{1}$ Cubic Double Perovskites, Phys. Rev. Lett. 118, 217202 (2017).
  * Natori  _et al._ (2018) W. M. H. Natori, E. C. Andrade, and R. G. Pereira, SU(4)-symmetric spin-orbital liquids on the hyperhoneycomb lattice, Phys. Rev. B 98, 195113 (2018).
  * Yamada  _et al._ (2018) M. G. Yamada, M. Oshikawa, and G. Jackeli, Emergent $\mathrm{SU}(4)$ Symmetry in $\alpha\text{$-$}{\mathrm{ZrCl}}_{3}$ and Crystalline Spin-Orbital Liquids, Phys. Rev. Lett. 121, 097201 (2018).
  * Ishikawa  _et al._ (2019) H. Ishikawa, T. Takayama, R. K. Kremer, J. Nuss, R. Dinnebier, K. Kitagawa, K. Ishii, and H. Takagi, Ordering of hidden multipoles in spin-orbit entangled $5{d}^{1}$ Ta chlorides, Phys. Rev. B 100, 045142 (2019).
  * de Farias  _et al._ (2020) C. S. de Farias, V. S. de Carvalho, E. Miranda, and R. G. Pereira, Quadrupolar spin liquid, octupolar Kondo coupling, and odd-frequency superconductivity in an exactly solvable model, Phys. Rev. B 102, 075110 (2020).
  * Chen and Wu (2021) G. Chen and C. Wu, Mott insulators with large local Hilbert spaces in quantum materials and ultracold atoms, arXiv:2112.02630 (2021).
  * Natori  _et al._ (2023) W. M. H. Natori, H.-K. Jin, and J. Knolle, Quantum liquids of the $S=3/2$ Kitaev honeycomb and related Kugel-Khomskii models, arXiv:2304.13378 (2023).
  * Li and Chen (2017) Y.-D. Li and G. Chen, Symmetry enriched U(1) topological orders for dipole-octupole doublets on a pyrochlore lattice, Phys. Rev. B 95, 041106 (2017).
  * Sibille  _et al._ (2020) R. Sibille, N. Gauthier, E. Lhotel, V. Porée, V. Pomjakushin, R. A. Ewings, T. G. Perring, J. Ollivier, A. Wildes, C. Ritter, T. C. Hansen, D. A. Keen, G. J. Nilsen, L. Keller, S. Petit, and T. Fennell, A quantum liquid of magnetic octupoles on the pyrochlore lattice, Nat. Phys. 16, 546 (2020).
  * Rayyan  _et al._ (2023) A. Rayyan, D. Churchill, and H.-Y. Kee, Field-induced Kitaev multipolar liquid in spin-orbit coupled ${d}^{2}$ honeycomb Mott insulators, Phys. Rev. B 107, L020408 (2023).
  * Khomskii and Streltsov (2021) D. I. Khomskii and S. V. Streltsov, Orbital effects in solids: Basics, recent progress, and opportunities, Chem. Rev. 121, 2992 (2021).
  * Streltsov and Khomskii (2020) S. V. Streltsov and D. I. Khomskii, Jahn-Teller Effect and Spin-Orbit Coupling: Friends or Foes? Phys. Rev. X 10, 031043 (2020).
  * Rau  _et al._ (2014) J. G. Rau, E. K.-H. Lee, and H.-Y. Kee, Generic Spin Model for the Honeycomb Iridates beyond the Kitaev Limit, Phys. Rev. Lett. 112, 077204 (2014).
  * Baker (1971) J. M. Baker, Interactions between ions with orbital angular momentum in insulators, Rep. Prog. Phys. 34, 109 (1971).
  * Gehring and Gehring (1975) G. A. Gehring and K. A. Gehring, Co-operative Jahn-Teller effects, Rep. Prog. Phys. 38, 1 (1975).
  * Voleti  _et al._ (2021) S. Voleti, A. Haldar, and A. Paramekanti, Octupolar order and Ising quantum criticality tuned by strain and dimensionality: Application to $d$-orbital Mott insulators, Phys. Rev. B 104, 174431 (2021).
  * Yadav  _et al._ (2018) R. Yadav, S. Rachel, L. Hozoi, J. van den Brink, and G. Jackeli, Strain- and pressure-tuned magnetic interactions in honeycomb Kitaev materials, Phys. Rev. B 98, 121107 (2018).
  * Santini  _et al._ (2009) P. Santini, S. Carretta, G. Amoretti, R. Caciuffo, N. Magnani, and G. H. Lander, Multipolar interactions in $f$-electron systems: The paradigm of actinide dioxides, Rev. Mod. Phys. 81, 807 (2009).
  * Iwahara  _et al._ (2018) N. Iwahara, V. Vieru, and L. F. Chibotaru, Spin-orbital-lattice entangled states in cubic ${d}^{1}$ double perovskites, Phys. Rev. B 98, 075138 (2018).
  * Wu (2008) C. Wu, Orbital Ordering and Frustration of $p$-Band Mott Insulators, Phys. Rev. Lett. 100, 200406 (2008).
  * Khaliullin  _et al._ (2021) G. Khaliullin, D. Churchill, P. P. Stavropoulos, and H.-Y. Kee, Exchange interactions, Jahn-Teller coupling, and multipole orders in pseudospin one-half $5{d}^{2}$ Mott insulators, Phys. Rev. Research 3, 033163 (2021).
  * Kubo and Kuramoto (2003) K. Kubo and Y. Kuramoto, Lattice Distortion and Octupole Ordering Model in CexLa1-xB6, J. Phys. Soc. Jpn. 72, 1859 (2003).
  * Kubo and Kuramoto (2004) K. Kubo and Y. Kuramoto, Octupole Ordering Model for the Phase IV of CexLa1-xB6, J. Phys. Soc. Jpn. 73, 216 (2004).
  * Paramekanti  _et al._ (2020) A. Paramekanti, D. D. Maharaj, and B. D. Gaulin, Octupolar order in $d$-orbital Mott insulators, Phys. Rev. B 101, 054439 (2020).
  * Churchill and Kee (2022) D. Churchill and H.-Y. Kee, Competing multipolar orders in a face-centered cubic lattice: Application to the osmium double perovskites, Phys. Rev. B 105, 014438 (2022).
  * Shiina  _et al._ (1997) R. Shiina, H. Shiba, and P. Thalmeier, Magnetic-Field Effects on Quadrupolar Ordering in a $\Gamma_{8}$-Quartet System CeB6, J. Phys. Soc. Jpn. 66, 1741 (1997).
  * Knolle  _et al._ (2014) J. Knolle, D. L. Kovrizhin, J. T. Chalker, and R. Moessner, Dynamics of a Two-Dimensional Quantum Spin Liquid: Signatures of Emergent Majorana Fermions and Fluxes, Phys. Rev. Lett. 112, 207203 (2014).
  * Fukui  _et al._ (2005) T. Fukui, Y. Hatsugai, and H. Suzuki, Chern Numbers in Discretized Brillouin Zone: Efficient Method of Computing (Spin) Hall Conductances, J. Phys Soc. Jpn. 74, 1674 (2005).
  * Zhang  _et al._ (2020) S.-S. Zhang, C. D. Batista, and G. B. Halász, Toward Kitaev’s sixteenfold way in a honeycomb lattice model, Phys. Rev. Res. 2, 023334 (2020).
  * Chulliparambil  _et al._ (2020) S. Chulliparambil, U. F. P. Seifert, M. Vojta, L. Janssen, and H.-H. Tu, Microscopic models for Kitaev’s sixteenfold way of anyon theories, Phys. Rev. B 102, 201111(R) (2020).
  * Fernandes  _et al._ (2019) R. M. Fernandes, P. P. Orth, and J. Schmalian, Intertwined Vestigial Order in Quantum Materials: Nematicity and Beyond, Annu. Rev. Condens. Matter Phys. 10, 133 (2019).
  * Murthy  _et al._ (1997) G. Murthy, D. Arovas, and A. Auerbach, Superfluids and supersolids on frustrated two-dimensional lattices, Phys. Rev. B 55, 3104 (1997).
  * Zhao and Liu (2008) E. Zhao and W. V. Liu, Orbital Order in Mott Insulators of Spinless $p$-Band Fermions, Phys. Rev. Lett. 100, 160403 (2008).
  * van den Brink  _et al._ (1998) J. van den Brink, W. Stekelenburg, D. I. Khomskii, G. A. Sawatzky, and K. I. Kugel, Elementary excitations in the coupled spin-orbital model, Phys. Rev. B 58, 10276 (1998).
  * Hou  _et al._ (2007) C.-Y. Hou, C. Chamon, and C. Mudry, Electron Fractionalization in Two-Dimensional Graphenelike Structures, Phys. Rev. Lett. 98, 186809 (2007).
  * A. C. Qu and P. Nigge and S. Link and G. Levy and M. Michiardi and P. L. Spandar and T. Matthé and M. Schneider and S. Zhdanovich and U. Starke and C. Gutiérrez and A. Damascelli (2022) A. C. Qu and P. Nigge and S. Link and G. Levy and M. Michiardi and P. L. Spandar and T. Matthé and M. Schneider and S. Zhdanovich and U. Starke and C. Gutiérrez and A. Damascelli, Ubiquitous defect-induced density wave instability in monolayer graphene, Sci. Adv. 8, eabm5180 (2022).
  * Jiang  _et al._ (2020) M.-H. Jiang, S. Liang, W. Chen, Y. Qi, J.-X. Li, and Q.-H. Wang, Tuning Topological Orders by a Conical Magnetic Field in the Kitaev Model, Phys. Rev. Lett. 125, 177203 (2020).
  * Zhang  _et al._ (2022) S.-S. Zhang, G. B. Halász, and C. D. Batista, Theory of the Kitaev model in a [111] magnetic field, Nat. Commun. 13, 399 (2022).
  * Lieb (1994) E. H. Lieb, Flux Phase of the Half-Filled Band, Phys. Rev. Lett. 73, 2158 (1994).