#  Eight-color chiral spin liquid in the $S=1$ bilinear-biquadratic model with
Kitaev interactions

Rico Pohle  Graduate School of Science and Technology, Keio University,
Yokohama 223-8522, Japan  Department of Applied Physics, University of Tokyo,
Hongo, Bunkyo-ku, Tokyo, 113-8656, Japan  Nic Shannon  Theory of Quantum
Matter Unit, Okinawa Institute of Science and Technology Graduate University,
Onna-son, Okinawa 904-0412, Japan  Yukitoshi Motome  Department of Applied
Physics, University of Tokyo, Hongo, Bunkyo-ku, Tokyo, 113-8656, Japan

(August 31, 2024)

###### Abstract

Multipolar spin systems provide a rich ground for the emergence of unexpected
states of matter due to their enlarged spin degree of freedom. In this study,
with a specific emphasis on $S=1$ magnets, we explore the interplay between
spin nematic states and spin liquids. Based on the foundations laid in the
prior work [R. Pohle et al., Phys. Rev. B 107, L140403 (2023)], we investigate
the $S=1$ Kitaev model with bilinear-biquadratic interactions, which
stabilizes, next to Kitaev spin liquid, spin nematic and triple-$q$ phases,
also an exotic chiral spin liquid. Through a systematic reduction of the spin
degree of freedom – from $\mathbb{CP}^{2}$ to $\mathbb{CP}^{1}$ and ultimately
to a discrete eight-color model – we provide an intuitive understanding of the
nature and origin of this chiral spin liquid. We find that the chiral spin
liquid is characterized by an extensive ground-state degeneracy, bound by a
residual entropy, extremely short-ranged correlations, a nonzero scalar spin
chirality marked by $\mathbb{Z}_{2}$ flux order, and a gapped continuum of
excitations. Our work contributes not only to the specific exploration of
$S=1$ Kitaev magnets but also to the broader understanding of the importance
of multipolar spin degree of freedom on the ground state and excitation
properties in quantum magnets.

##  I Introduction

Multipolar systems exhibit higher-order moments, such as quadrupoles and
octupoles, even in absence of monopoles and dipoles. The inherent fluctuations
and allowed higher-order interactions in these systems provide a rich
landscape of unconventional phases, evident from multipolar order observed in
quantum magnets [1, 2, 3, 4, 5], unconventional superconductors [6, 7, 8, 9,
10, 11], and topological insulators [12, 13]. Particular exotic examples
include multipolar quantum spin ice [14, 15, 16], $\mathbb{CP}^{2}$ skyrmion
crystals [17, 18, 19], and unconventional orders in the higher-spin Kondo
lattice model [20, 21], demonstrating the diverse range of intriguing quantum
phenomena present in multipolar condensed matter systems.

Multipolar magnets made of $S=1$ moments are particularly intriguing, as they
strike a delicate balance with a spin length that is “small enough” to exhibit
strong quantum effects, while simultaneously being “large enough” to give rise
to onsite quadrupole moments [22, 23, 24, 25]. Quadrupoles break spin-rotation
symmetry by selecting an axis without specifying a particular direction,
closely resembling properties observed in classical liquid crystals [26, 27,
28]. This unique characteristic makes them an ideal playground for exploring
topological defects [29, 30], investigating out-of-equilibrium effects [31,
32], and drawing potential analogies to gravity [33, 34, 35]. Despite being
theoretically proposed almost 60 years ago [36, 37], the experimental
exploration of $S=1$ magnets has faced substantial challenges, primarily since
their quadrupolar order remains hidden from conventional magnetic probes like
neutron scattering [38, 39]. Fortunately, recent advancements in experimental
techniques such as Raman scattering [40, 41], resonant inelastic X-ray
scattering (RIXS) [42, 43, 44], and the exploration of physical properties
like the magnetocaloric effect [45], combined with the development of powerful
numerical techniques [46, 47, 48, 49, 50, 51], have expanded our capabilities
to interpret thermodynamic and dynamic properties of multipolar magnets in
real materials. The successful application of these diverse methods has
significantly advanced our understanding of various $S=1$ compounds, including
NiGa2S4 [52, 53, 41], FeI2 [54, 55], and Ba2FeSi2O7 [56].

Recently, higher-spin Kitaev systems have attracted substantial attention [57,
58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68], as they offer a new avenue for
studying quantum spin liquids with potentially multipolar character [69, 70]
and fractionalized excitations carrying multipolar moments [66]. Theoretical
proposals suggest that magnets with strong Hund’s coupling and spin-orbit
interactions could pave the way for experimental realizations [71], with
promising candidate materials such as NaNi2BiO6-δ [72] and KNiAsO4 [73].
Furthermore, $S=1$ magnets naturally allow for higher-order interactions, such
as biquadratic interactions, which effectively couple the quadrupole degree of
freedom between individual spins. In this context, the recent work by the
authors explored the intricate interplay between spin nematics and spin
liquids in the $S=1$ Kitaev model with bilinear-biquadratic (BBQ) interactions
[74]. The richness of this model unfolds with a plethora of phases, including
spin nematic, and Kitaev spin liquid states, alongside triple-$q$ order with
nonzero scalar spin chirality, ferro, antiferro, zigzag and stripy phases.
Particularly intriguing is the stabilization of exotic phases where dipole and
quadrupole moments compete. These include a quasi-one-dimensional coplanar
phase, a twisted conical phase, and a noncoplanar ordered state. A
comprehensive phase diagram showcasing these diverse phases is illustrated in
Fig. 1(a). Interestingly, in the vicinity of the noncoplanar ordered state,
the model hosts an unconventional chiral spin liquid (CSL) at finite
temperature, which exhibits an extensive degeneracy of states, retains the
Kitaev spin liquid (SL) feature of $\mathbb{Z}_{2}$ flux order and shows a
nonzero scalar spin chirality [74].

Figure 1:  Ground-state phase diagrams in the 2D equirectangular projection as
function of $\theta$ and $\phi$ [see Eq. (2)]. We plot the second derivative
of the internal energy as contour color to make phase boundaries visible. (a)
Solution of $\mathcal{H}^{\mathcal{A}}_{{\sf BBQ-K}}$ [see Eq. (7)] in spin
space $\mathbb{CP}^{2}$, which accounts for all local degrees of freedom of an
$S=1$ moment, including dipolar and quadrupolar components. (b) Solution of
$\mathcal{H}^{\mathcal{S}}_{{\sf BBQ-K}}$ [see Eq. (1)] in the spin space
$\mathbb{CP}^{1}$, which considers only dipolar degree of freedom of classical
Heisenberg spins. (c) Solution of $\mathcal{H}^{{\sf 8c}}_{{\sf BBQ-K}}$ for
the discrete eight-color model [see Eq. (15)], where the dipolar degree of
freedom are discretized and allowed to point only at the corners of a unit
cube (see definition of spin states in Table 1). (d) Real-space configurations
of selected phases. Ferromagnetic (FM), antiferromagnetic (AFM), and
triple-$q$ chiral magnetic orders (FM 3$Q$ chiral and AFM 3$Q$ chiral), as
well as the semiclassical analog of the Kitaev spin liquid (SL), are
stabilized in all three phase diagrams. Noncoplanar (NC), (quasi-)one-
dimensional (q1D / 1D), canted planar (CPL), and twisted conical (TC) phases
are stabilized for the models with continuous degrees of freedom in (a) and
(b), while ferroquadrupolar (FQ), and semiordered (SO) phases are purely
quadrupolar in nature and are only stabilized in (a). FM and AFM chiral spin
liquid (CSL) states emerge as the ground state along two singular lines in (b)
and form a phase in (c), while they are not stabilized in the ground state in
(a) but become stable at finite temperatures (see discussion in Sec. III.4).
Results in (a) and (b) were obtained by variational energy minimization for a
finite-size cluster of $N_{\rm S}=28\ 800$ spins under periodic boundary
conditions, while (c) was solved analytically for the thermodynamic limit.

In this article, we expand on the theory discussed in our prior work, Ref.
[74], with the goal to present a pedagogical explanation of the inherent
nature of the CSL. To achieve this, we start by representing local $S=1$
magnetic moments as SU(3) spin-coherent states in the spin space
$\mathbb{CP}^{2}$. In this representation, spins exhibit both onsite
quadrupole and dipole characteristics, as expected for $S=1$ moments [47]. We
proceed by systematically reducing the degrees of freedom for local spins to
SU(2) spin-coherent states in the spin space $\mathbb{CP}^{1}$. This
representation is characterized by classical Heisenberg spins that exclusively
contain dipole moments. Finally, we arrive at a minimal model, the eight-color
model, where dipolar spins are confined to eight discrete spin states. The
corresponding phase diagrams for each of these cases are depicted in Figs.
1(a), 1(b), and 1(c). Many of the most interesting properties of the CSL are
captured by the eight-colour model, which allows us to understand its physical
properties analytically. We clarify the origin of the bond constraints that
unequivocally determine the nature of this CSL, and validate our findings
through classical Monte Carlo simulations. Furthermore, we present the
dynamical properties of this phase, unveiling a gapped continuum of
excitations.

The remainder of this article is structured as follows. In Sec. II, we
introduce the $S=1$ BBQ-Kitaev model on the honeycomb lattice. We analyze the
same model for different local spin degrees of freedom, namely, SU(3) spin-
coherent states of $\mathbb{CP}^{2}$ (the $\mathbb{CP}^{2}$ model) which
contain onsite dipole and quadrupole moments, SU(2) spin-coherent states of
$\mathbb{CP}^{1}$ (the $\mathbb{CP}^{1}$ model) which contain only onsite
dipole moments, and the eight-color model which allows for only eight
discretized dipolar spin states. In Sec. III, we present ground-state and
finite-temperature results of the BBQ-Kitaev model. In Sec. III.1, we discuss
the influence of the local spin degree of freedom by explicitly comparing the
ground-state phase diagrams between the $\mathbb{CP}^{2}$, the
$\mathbb{CP}^{1}$, and the eight-color models. In Sec. III.2, we discuss the
eight-color model and provide an analytical explanation for all fundamental
properties of the observed CSL, substantiated by comparisons to numerical
simulations. In Sec. III.3, we show numerical results for the
$\mathbb{CP}^{1}$ model, and show that thermodynamic properties of the eight-
color CSL remain intact even for continuous degrees of freedom. In Sec. III.4,
we show numerical results for the $\mathbb{CP}^{2}$ model, demonstrating that
the CSL becomes an entropy-driven spin liquid at finite temperature. In Sec.
III.5, we present dynamical signatures of the CSL in the $\mathbb{CP}^{2}$
model, showing a largely gapped continuum of excitations. We conclude our work
in Sec. IV with a brief summary and perspectives for possible future
directions and experimental realizations.

##  II BBQ-Kitaev model with different degrees of freedom

We are interested in solving the bilinear-biquadratic (BBQ) model under
influence of Kitaev interactions for $S=1$ magnetic moments on the honeycomb
lattice. The Hamiltonian is given by

$\displaystyle\mathcal{H}^{\mathcal{S}}_{{\sf BBQ-K}}=\sum_{\langle
i,j\rangle}$ $\displaystyle\left[\ J_{1}\ {\bf S}_{i}\cdot{\bf S}_{j}\ +\
J_{2}\ ({\bf S}_{i}\cdot{\bf S}_{j})^{2}\ \right]$ (1) $\displaystyle\
+K\sum_{\alpha=x,y,z}\ \sum_{\langle
ij\rangle_{\alpha}}S_{i}^{\alpha}S_{j}^{\alpha}\,,$

where $J_{1}$, $J_{2}$, and $K$ respectively account for the Heisenberg
(bilinear), biquadratic, and Kitaev interaction strengths on nearest-neighbor
bonds. The index $\alpha=x,y,z$ selects the spin- and bond-anisotropic Kitaev
interactions on the honeycomb lattice [75]. We normalize the total interaction
strength as

$\left(J_{1},J_{2},K\right)=\left(\sin{\theta}\cos{\phi},\sin{\theta}\sin{\phi},\cos{\theta}\right)\,,$
(2)

and use the angles $\theta$ and $\phi$ as new model parameters. In this form,
$\mathcal{H}^{\mathcal{S}}_{{\sf BBQ-K}}$ recovers the well-known limits of
the BBQ model for $\theta/\pi=0.5$ ($K=0$), the AFM Kitaev model at
$\theta/\pi=0$ ($K=1$, $J_{1}=J_{2}=0$), the FM Kitaev model at $\theta/\pi=1$
($K=-1$, $J_{1}=J_{2}=0$), and the Kitaev-Heisenberg model for $\phi/\pi=0$
and $1$ ($J_{2}=0$).

An $S=1$ moment contains three magnetic states with $S^{z}=\\{1,0,-1\\}$,
which are mathematically described by the $\mathfrak{su}(3)$ algebra and
contain magnetic dipole and quadrupole components. To correctly describe all
fluctuations of such an $S=1$ moment, we adopt the recently developed U(3)
formalism [47], an approach equivalent to the concept of SU(3) spin-coherent
states [22, 76, 77, 78, 46, 48, 49], and suitable to access not only
thermodynamic properties, but also dynamical properties of $S=1$ magnets at
finite temperatures. By embedding the underlying $\mathfrak{su}(3)$ algebra
into the larger $\mathfrak{u}(3)$ algebra with an additional spin-length
constraint, this formalism becomes capable of simultaneously evaluating dipole
and quadrupole fluctuations at each site, despite the drawback of losing
quantum entanglement across the lattice.

In the U(3) formalism, we express each localized $S=1$ moment by a $3\times 3$
Hermitian matrix

$\mathcal{A}^{\alpha}_{\beta}=({\rm d}^{\alpha})^{*}\ {\rm d}_{\beta}\,,$ (3)

which we shall call henceforth $\mathcal{A}$ matrix. These $\mathcal{A}$
matrices can be written in terms of complex vectors ${\bf d}$

${\bf d}=\begin{pmatrix}x_{1}+i\ x_{2}\\\ x_{3}+i\ x_{4}\\\ x_{5}+i\
x_{6}\end{pmatrix},$ (4)

also known as directors [24, 23]. We implicitly respect the spin-length
constraint on the directors

${\bf d}^{*}{\bf d}=|{\bf d}|^{2}=1\,,$ (5)

by parametrizing its components as

$\displaystyle x_{1}$ $\displaystyle=\theta_{2}^{1/4}\ \theta_{1}^{1/2}\
\sin{\phi_{1}}\ ,$ (6) $\displaystyle x_{2}$ $\displaystyle=\theta_{2}^{1/4}\
\theta_{1}^{1/2}\ \cos{\phi_{1}}\ ,$ $\displaystyle x_{3}$
$\displaystyle=\theta_{2}^{1/4}\ (1-\theta_{1})^{1/2}\ \sin{\phi_{2}}\ ,$
$\displaystyle x_{4}$ $\displaystyle=\theta_{2}^{1/4}\ (1-\theta_{1})^{1/2}\
\cos{\phi_{2}}\ ,$ $\displaystyle x_{5}$
$\displaystyle=(1-\theta_{2}^{1/2})^{1/2}\ \sin{\phi_{3}}\ ,$ $\displaystyle
x_{6}$ $\displaystyle=(1-\theta_{2}^{1/2})^{1/2}\ \cos{\phi_{3}}\ ,$

with $0\leq\theta_{1},\theta_{2}\leq 1$ and
$0\leq\phi_{1},\phi_{2},\phi_{3}<2\pi$. To fix the unphysical gauge we choose
to make the $z$ component of ${\bf d}$ purely real by setting
$\phi_{3}=\pi/2$, providing in total four degrees of freedom with a restrained
order-parameter space, formally known as the complex projective plane
$\mathbb{CP}^{2}$. In this formalism, Eq. (1) can be rewritten in a bilinear
form as

$\displaystyle\mathcal{H}^{\mathcal{A}}_{{\sf BBQ-K}}$ $\displaystyle=$ (7)
$\displaystyle\sum_{\langle ij\rangle}\Big{[}J_{1}$
$\displaystyle\mathcal{A}^{\alpha}_{i\beta}\mathcal{A}^{\beta}_{j\alpha}+(J_{2}-J_{1})\mathcal{A}^{\alpha}_{i\beta}\mathcal{A}^{\alpha}_{j\beta}+J_{2}\mathcal{A}^{\alpha}_{i\alpha}\mathcal{A}^{\beta}_{j\beta}\Big{]}$
$\displaystyle-K\sum_{\langle ij\rangle_{\alpha}}\epsilon^{\alpha\ \gamma}_{\
\beta}\epsilon^{\alpha\ \eta}_{\
\delta}\mathcal{A}^{\beta}_{i\gamma}\mathcal{A}^{\delta}_{j\eta}\,,$

where $\epsilon^{\alpha\ \gamma}_{\ \beta}$ is the Levi-Civita symbol; here
and hereafter, we adopt the Einstein convention of summing over repeated
indices. We shall refer to this form of the Hamiltonian in the following as
“$\mathbb{CP}^{2}$ model”. The matrix $\mathcal{A}_{i}$ simultaneously
incorporates the information of dipole and quadrupole moments at site $i$,
which can be extracted respectively with

$\displaystyle S_{i}^{\alpha}$ $\displaystyle=-i\epsilon^{\alpha\ \gamma}_{\
\beta}\mathcal{A}^{\beta}_{i\gamma}\,,$ (8) $\displaystyle
Q^{\alpha\beta}_{i}$
$\displaystyle=-\mathcal{A}^{\alpha}_{i\beta}-\mathcal{A}^{\beta}_{i\alpha}+\frac{2}{3}\delta^{\alpha\beta}\mathcal{A}^{\gamma}_{i\gamma}\,,$
(9)

where $\delta^{\alpha\beta}$ is the Kronecker delta. Here, we express spin
quadrupoles in their symmetric and traceless rank-2 tensor form.

To compute the averaged spin-dipole and spin-quadrupole norms we respectively
use

$\displaystyle|{\bf S}|$ $\displaystyle=\frac{1}{N_{\rm S}}\sum^{N_{\rm
S}}_{i}|{\bf S}_{i}|\,,$ (10) $\displaystyle|{\bf Q}|$
$\displaystyle=\frac{1}{N_{\rm S}}\sum^{N_{\rm S}}_{i}|{\bf Q}_{i}|\,,$ (11)

where vector components are expressed as

${\bf S}_{i}=\begin{pmatrix}S_{i}^{x}\\\ S_{i}^{y}\\\
S_{i}^{z}\end{pmatrix}\,,$ (12)

and

${\bf Q}_{i}=\begin{pmatrix}Q_{i}^{x^{2}-y^{2}}\\\ Q_{i}^{3z^{2}-r^{2}}\\\
Q_{i}^{xy}\\\ Q_{i}^{xz}\\\
Q_{i}^{yz}\end{pmatrix}=\begin{pmatrix}\frac{1}{2}\left(Q_{i}^{xx}-Q_{i}^{yy}\right)\\\
\frac{1}{\sqrt{3}}\left(Q_{i}^{zz}-\frac{1}{2}\left(Q_{i}^{xx}+Q_{i}^{yy}\right)\right)\\\
Q_{i}^{xy}\\\ Q_{i}^{xz}\\\ Q_{i}^{yz}\\\ \end{pmatrix}\,.$ (13)

In Sec. III.1, we find many nontrivial phases for the $\mathbb{CP}^{2}$ model
in Eq. (7), especially when the interactions frustrate the model. To gain a
better intuition of the underlying physics in such nontrivial phases, we
investigate Eq. (1) also for the spin space $\mathbb{CP}^{1}$, by effectively
excluding all quadrupole fluctuations and allowing only for fluctuations of
local dipole moments. This simplification corresponds to the treatment of
classical Heisenberg spins, formally expressed as spin-coherent states of
SU(2), which are parametrized by a point on the Bloch sphere

${\bf S}=\begin{pmatrix}\sin{\theta_{1}}\cos{\phi_{1}}\\\
\sin{\theta_{1}}\sin{\phi_{1}}\\\ \cos{\theta_{1}}\end{pmatrix},$ (14)

with $0\leq\theta_{1}\leq 1$ and $0\leq\phi_{1}<2\pi$, the two local degrees
of freedom for each spin on the lattice. In what follows, we refer to this
model as “$\mathbb{CP}^{1}$ model”.

| bright | dark
---|---|---
green | $\boldsymbol{\sigma}_{\rm g}=\frac{1}{\sqrt{3}}\left\\{-1,-1,-1\right\\}$ | $\bar{\boldsymbol{\sigma}}_{\rm g}=\frac{1}{\sqrt{3}}\left\\{+1,+1,+1\right\\}$
yellow | $\boldsymbol{\sigma}_{\rm y}=\frac{1}{\sqrt{3}}\left\\{+1,+1,-1\right\\}$ | $\bar{\boldsymbol{\sigma}}_{\rm y}=\frac{1}{\sqrt{3}}\left\\{-1,-1,+1\right\\}$
blue | $\boldsymbol{\sigma}_{\rm b}=\frac{1}{\sqrt{3}}\left\\{-1,+1,+1\right\\}$ | $\bar{\boldsymbol{\sigma}}_{\rm b}=\frac{1}{\sqrt{3}}\left\\{+1,-1,-1\right\\}$
red | $\boldsymbol{\sigma}_{\rm r}=\frac{1}{\sqrt{3}}\left\\{+1,-1,+1\right\\}$ | $\bar{\boldsymbol{\sigma}}_{\rm r}=\frac{1}{\sqrt{3}}\left\\{-1,+1,-1\right\\}$

Table 1:  Definition of spin directions in the eight-color model of Eq. (15).
For pedagogical reasons, we assign four colors to spins, namely “green”,
“yellow”, “blue”, and “red”, and distinguish between their bright and dark
contrast. Directions are explicitly shown in (a) for bright spins, and in (b)
for dark spins.

Moreover, to transparently understand the nature of exotic chiral ordered and
disordered phases (see Fig. 1 and detailed discussion in Sec. III.2), we even
further simplify Eq. (1) into the “eight-color model”, whose Hamiltonian is
given by

$\displaystyle\mathcal{H}^{{\sf 8c}}_{{\sf BBQ-K}}=\sum_{\langle i,j\rangle}$
$\displaystyle\left[\ J_{1}\
\boldsymbol{\sigma}_{i}\cdot\boldsymbol{\sigma}_{j}\ +\ J_{2}\
(\boldsymbol{\sigma}_{i}\cdot\boldsymbol{\sigma}_{j})^{2}\ \right]$ (15)
$\displaystyle\ +K\sum_{\alpha=x,y,z}\ \sum_{\langle
ij\rangle_{\alpha}}\sigma_{i}^{\alpha}\sigma_{j}^{\alpha}\,,$

where discrete, noncoplanar dipole spins $\boldsymbol{\sigma}$ are only
allowed to point along the eight corners of a unit cube. As shown in Table 1,
we assign colors to spins such as “green”, “yellow”, “blue”, and “red”, and
distinguish between “bright” spins, for $\boldsymbol{\sigma}_{\rm g}$,
$\boldsymbol{\sigma}_{\rm y}$, $\boldsymbol{\sigma}_{\rm b}$, and
$\boldsymbol{\sigma}_{\rm r}$, and “dark” spins, for
$\bar{\boldsymbol{\sigma}}_{\rm g}$, $\bar{\boldsymbol{\sigma}}_{\rm y}$,
$\bar{\boldsymbol{\sigma}}_{\rm b}$, and $\bar{\boldsymbol{\sigma}}_{\rm r}$.

##  III Results

In our previous work in Ref. [74], we studied the semiclassical ground state
and thermodynamic properties of $\mathcal{H}^{\mathcal{S}}_{{\sf BBQ-K}}$ in
Eq. (1), and revealed the existence of a chiral spin liquid (CSL) at finite
temperatures. To build up an intuition for the microscopic origin and nature
of this exotic CSL state, we systematically restrict the local spin degree of
freedom in the present study. We begin with the $\mathfrak{su}(3)$
representation of Eq. (1), the $\mathbb{CP}^{2}$ model, given by
$\mathcal{H}^{\mathcal{A}}_{{\sf BBQ-K}}$ in Eq. (7), which captures all
allowed fluctuations for an $S=1$ moment, namely continuous dipolar and
quadrupolar fluctuations described in the spin space $\mathbb{CP}^{2}$.
Further, we solve the $\mathbb{CP}^{1}$ model, namely the original model
$\mathcal{H}^{\mathcal{S}}_{{\sf BBQ-K}}$ in Eq. (1) only with continuous
dipolar fluctuations described by the spin space $\mathbb{CP}^{1}$ [see Eq.
(14)]. Finally, we consider the discrete eight-color model given by
$\mathcal{H}^{{\sf 8c}}_{{\sf BBQ-K}}$ in Eq. (15), where dipolar spin moments
are allowed to point only along eight discrete directions (see Table 1). All
three models describe the same BBQ-Kitaev interactions, however for different
degrees of freedom in spin space. While certain results have been previously
presented in our earlier work [74], we will include some of them here to make
the paper self-contained.

###  III.1 Ground-state phase diagrams

In Figs. 1(a)–1(c) we respectively show the ground-state phase diagrams of
$\mathcal{H}^{\mathcal{A}}_{{\sf BBQ-K}}$, $\mathcal{H}^{\mathcal{S}}_{{\sf
BBQ-K}}$, and $\mathcal{H}^{{\sf 8c}}_{{\sf BBQ-K}}$, solved in the spin space
$\mathbb{CP}^{2}$, $\mathbb{CP}^{1}$, and for discrete eight-color spins. To
obtain Figs. 1(a) and 1(b), we performed variational energy minimization,
using the gradient descent method based on the machine learning library JAX
[79, 80] for sufficiently large finite-size clusters of $N_{\rm S}=2L^{2}=28\
800$ spins (linear dimension $L=120$) under periodic boundary conditions (see
Appendix A for further details). Figure 1(c) was obtained in the thermodynamic
limit from analytic comparison of local bond energies. We plot the phase
diagrams in their 2D equirectangular projections parametrized by Eq. (2), and
visualize phase boundaries by plotting
$\partial^{2}E/\partial\theta^{2}+\partial^{2}E/\partial\phi^{2}$, the sum of
the second-order derivatives of the internal energy $E$ with respect to
$\theta$ and $\phi$. The corresponding real-space spin configurations of
selected states are shown in Fig. 1(d) with dipole moments depicted as arrows,
and quadrupole moments in form of “directors” with “doughnut-shaped” spin-
probability distributions [24, 23]. It is worth noting that the phase diagrams
for $\mathbb{CP}^{1}$ and eight-color models exhibit self-duality around the
high-symmetry points $(\phi/\pi,\theta/\pi)=(0.5,0.5)$ and
$(\phi/\pi,\theta/\pi)=(1.5,0.5)$. In contrast, the phase diagram for
$\mathbb{CP}^{2}$ lacks this self-dual symmetry, due to the presence of
quadrupole moments.

####  III.1.1 Phase diagram for $\mathbb{CP}^{2}$ model

In Fig. 1(a) we show the ground-state phase diagram of
$\mathcal{H}^{\mathcal{A}}_{{\sf BBQ-K}}$ in Eq. (7), by considering the full
spin order-parameter space $\mathbb{CP}^{2}$. When the biquadratic interaction
is ferromagnetic, namely $J_{2}$ is negative $(1.0<\phi/\pi<2.0)$, the model
offers a combination of phases which have previously been reported on the
$S=1$ Kitaev-Heisenberg model [71, 61] and the BBQ model [81]. These comprise
dipolar orders, in the form of ferromagnetic (FM), antiferromagnetic (AFM),
zigzag, and stripy phases, together with a ferroquadrupolar (FQ) spin nematic
phase.

In contrast, when $J_{2}$ is positive $(0.0<\phi/\pi<1.0)$, the model provides
a number of unconventional phases which emerge from the competition between
frustrated interactions. Adjacent to the zigzag and stripy phases, we find
noncoplanar chiral ordered phases, also known as “tetrahedral” and “cubic”
states [82, 83, 84], in which spins point to the corners of a unit cube [see
Fig. 1(d)]. These chiral configurations are represented by superpositions of
three spiral states, and hence called triple-$q$ states. The FM 3$Q$ chiral
state is characterized by an extended eight-site magnetic unit cell, where
four out of eight possible spin directions are selected, showing a net scalar
spin chirality of $|\kappa|^{\rm FM\ chiral}=8/(3\sqrt{3})$ [see the
definition of $\kappa$ in Eq. (60) of Appendix B]. The AFM 3$Q$ chiral state
is essentially a reflection of the FM 3$Q$ chiral state, achieved by inverting
spins on one of the two sublattices on the honeycomb lattice [see Fig. 1(d)],
giving $|\kappa|^{\rm AFM\ chiral}=16/(3\sqrt{3})$. We note that the stripy
and FM 3$Q$ chiral (zigzag and AFM 3$Q$ chiral) states are energetically
degenerate at $\phi/\pi=0$ ($\phi/\pi=1$) along the vertical line for
$0.75\leq\theta/\pi<1.0$ ($0.0<\theta/\pi\leq 0.25$). However, this degeneracy
is lifted by thermal fluctuations, which select the stripy (zigzag) phase, as
confirmed by our MC simulations (not shown).

Between FM, AFM, and 3$Q$ chiral ordered phases, dipole and quadrupole moments
mix in nontrivial ways, and produce a range of exotic states. In the limit of
the BBQ model ($\theta/\pi=0.5$), our semiclassical method realizes a
semiordered (SO) (or semidisordered) [22, 85], purely quadrupolar state,
connecting the two SU(3) points at $\phi/\pi=0.25$ and $\phi/\pi=0.5$ 111 We
notice several complicated phases around the two SU(3) points at
$(\theta/\pi,\phi/\pi)=(0.5,0.25)$ and $(0.5,0.5)$, which will be discussed
elsewhere. . We note that this state is expected to be replaced by a plaquette
valence bond crystal when quantum entanglement is fully taken into account
[87]. Introducing nonzero Kitaev interactions ($\theta/\pi\neq 0.5$)
immediately induces the formation of a nonzero dipole moment $|{\bf S}|\neq
0$. To minimize the dominant BBQ interactions, spins prefer to align
orthogonal to their nearest neighbors within the $xy$, $yz$, or $zx$ plane. In
such a configuration, the Kitaev interaction energy is minimized on two bonds,
while leaving the third one to be zero. Hence, dominant correlations prevail
along zigzag chains and stabilize coplanar, quasi-one-dimensional (q1D) states
in extended regions away from the BBQ limit. However, the presence of small
quadrupolar correlations will induce a weak interchain coupling to form a two-
dimensional order.

Sandwiched between the q1D coplanar phases and the AFM/FM ordered phases,
$\mathcal{H}^{\mathcal{A}}_{{\sf BBQ-K}}$ minimizes its energy by forming spin
textures with noncoplanar (NC) orientation and incommensurate twisted conical
(TC) states. Theses ordered states will give way to a finite-temperature
chiral spin liquid (CSL) – not shown in Fig. 1(a) – which shall be the major
focus of this study. We will discuss signatures of the CSL state explicitly
for $\mathcal{H}^{\mathcal{A}}_{{\sf BBQ-K}}$ by including all allowed spin
fluctuations of $\mathbb{CP}^{2}$ in Sec. III.4, while building up an
intuitive understanding of its microscopic origin by systematically
restricting the allowed spin fluctuations to $\mathbb{CP}^{1}$ and further to
discretized states only in Secs. III.3 and III.2, respectively.

####  III.1.2 Phase diagram for $\mathbb{CP}^{1}$ model

Many phases in Fig. 1(a) show dominant dipolar order mainly coming from the
presence of bilinear Kitaev and Heisenberg interactions. To build up an
intuitive understanding of the microscopic origin of these phases, and to
quantitatively understand the influence of local quantum fluctuations in our
model, we are going to restrict the order parameter space from
$\mathbb{CP}^{2}$ to $\mathbb{CP}^{1}$, by expressing an $S=1$ moment by a
classical Heisenberg vector, as parametrized by Eq. (14). In this way we only
allow for dipolar spin fluctuations and exclude onsite quadrupole moments.

In this framework, by minimizing the energy of
$\mathcal{H}^{\mathcal{S}}_{{\sf BBQ-K}}$ in Eq. (1), we obtain the ground-
state phase diagram in Fig. 1(b). For negative $J_{2}$ $(1.0<\phi/\pi<2.0)$
the model captures the same dominant dipolar FM, AFM, zigzag, and stripy
ordered phases, as in the $\mathbb{CP}^{2}$ case. However, due to the
exclusion of onsite quadrupolar degree of freedom, the FQ state is absent and
has been replaced by a singular point at $\phi/\pi=1.5$ and $\theta/\pi=0.5$
where FM, AFM, zigzag, and stripy phases are energetically degenerate.

In the region where $J_{2}$ is positive $(0.0<\phi/\pi<1.0)$, the model
reproduces many phases as seen in Fig. 1(a), while the phase boundaries are
rather strongly modified. Also here, we find noncoplanar FM and AFM 3$Q$
chiral phases with same physical properties similar to the $\mathbb{CP}^{2}$
case. Between AFM/FM ordered phases and AFM/FM 3$Q$ chiral phases the model
stabilizes NC phases and canted planar (CPL) phases with spin configurations
shown in Fig. 1(d). Here, the NC phases are purely made of spin dipoles and
exist in much wider regions compared to the $\mathbb{CP}^{2}$ case. At the
phase boundary between the NC and CPL phases ($\phi/\pi=0.5$), the model
stabilizes coplanar, truly 1D phases, made of decoupled zigzag chains with a
subextensive ground-state manifold [see Fig. 1(d)]. While these 1D phases
appear on a singular line in the $\mathbb{CP}^{1}$ model, the presence of
quadrupole fluctuations in the $\mathbb{CP}^{2}$ case stabilizes q1D phases
(not truly 1D due to presence of small but nonzero quadrupolar order), in much
wider regions of the phase diagram [see Fig. 1(a)].

Furthermore, the absence of quadrupole degree of freedom allows the CSL state
to exist in the ground state at the singular lines

$\frac{\phi}{\pi}=\frac{1}{\pi}\arctan\left(\frac{3}{2}\right)\approx
0.312833\,,\quad 0.0<\frac{\theta}{\pi}<0.5\,,$ (16)

for the FM CSL and at

$\frac{\phi}{\pi}=1-\frac{1}{\pi}\arctan\left(\frac{3}{2}\right)\approx
0.687167\,,\quad 0.5<\frac{\theta}{\pi}<1.0\,,$ (17)

for the AFM CSL, deep inside the NC phase. The analytic derivation of these
values is described in Appendix C.

In the CSL ground state, spin-dipole moments point exactly to the eight
discrete corners of the unit cube, similar to the ordered AFM/FM 3$Q$ chiral
states [see Fig. 1(d)]. However, in strong contrast to the chiral ordered
states, the CSL does not break translational symmetries of the lattice, and
hence exhibits an extensive degeneracy of states. This motivates us to further
simplify the model by restricting the continuous spin space of
$\mathbb{CP}^{1}$ to only eight discrete states (eight-color model), where
dipolar spins are allowed to point only along discretized corners of the unit
cube. We will explore this model in detail in the following subsections, and
discuss ground-state and thermodynamic properties of the CSL explicitly for
$\mathcal{H}^{\mathcal{S}}_{{\sf BBQ-K}}$ on the $\mathbb{CP}^{1}$ model in
Sec. III.3.

####  III.1.3 Phase diagram for eight-color model

To understand the intriguing nature of this CSL, we analyze the eight-color
model described by $\mathcal{H}^{{\sf 8c}}_{{\sf BBQ-K}}$ in Eq. (15), which
enables us to extract its most essential properties analytically. In this
model we restrict the spin parameter space by allowing for only eight discrete
dipolar spin directions, pointing along the corners of a unit cube, as defined
in Table 1.

As shown in the ground-state phase diagram in Fig. 1(c), this model stabilizes
trivial magnetically-ordered states such as FM and AFM states in the region
where the biquadratic interaction $J_{2}$ is mostly negative
$(1.0<\phi/\pi<2.0)$. For predominantly positive $J_{2}$ $(0.0<\phi/\pi<1.0)$,
the model stabilizes nontrivial chiral ordered, triple-$q$ states, as also
observed in the $\mathbb{CP}^{2}$ and $\mathbb{CP}^{1}$ models in Figs. 1(a)
and 1(b), respectively. Since coplanar spin configurations are not allowed in
this model, zigzag, stripy, 1D, and other states with canting angles away from
the (111) axes are absent.

Importantly, between trivial magnetically-ordered and triple-$q$ states, the
model stabilizes the CSL states in wide regions of the phase diagram. The FM
CSL phase exists for

$\frac{1}{\pi}\ \arctan{\left(\frac{3}{4}\right)}<\frac{\phi}{\pi}<0.5\,,\quad
0.0<\frac{\theta}{\pi}<0.5\,,$ (18)

while the AFM CSL phase appears for

$0.5<\frac{\phi}{\pi}<1-\frac{1}{\pi}\arctan{\left(\frac{3}{4}\right)}\,,\quad
0.5<\frac{\theta}{\pi}<1.0\,.$ (19)

The phase boundaries can be obtained analytically by simple energy comparison,
as described in Appendix C.

Both states are classical spin liquids with a nonzero scalar spin chirality
and residual entropy at zero temperature, which is associated with their
extensive number of degenerate ground states. We note that the phase boundary
lines between chiral order and disordered phases at $\phi/\pi=0.5$ and
$\theta/\pi=0.5$ also host classical spin liquids, which, however, has
different properties than the CSL discussed here.

By comparison between measured physical observables, we argue that this CSL
with all its properties survives in the $\mathbb{CP}^{1}$ and
$\mathbb{CP}^{2}$ models, while additional degrees of freedom in spin space
only weaken its stability in the phase diagram. For pedagogical reasons, we
will start our journey by investigating the CSL on the eight-color model.

###  III.2 Analysis of eight-color model

In Fig. 1(c), we see that the eight-color model stabilizes the CSL ground
state in wide regions of the phase diagram [see Eqs. (18) and (19)]. In Sec.
III.2.1, we use the simplicity of the eight-color model to analytically obtain
bond-dependent spin constraints which determine the ground-state properties of
this CSL. In Sec. III.2.2, we validate our analytical results by numerical MC
simulations and compare thermodynamic properties between the CSL and chiral
ordered states.

####  III.2.1 Eight-color chiral spin liquid

The eight-color model given by $\mathcal{H}^{{\sf 8c}}_{{\sf BBQ-K}}$ in Eq.
(15) allows us to understand the nature of the CSL analytically from simple
comparison of local bond energies. In the following, we consider the FM CSL
for parameters outlined in Eq. (18), where all couplings are positive, with
$0<J_{1}<0.8$, $J_{2}>0.6$, and $K>0$. Since $\mathcal{H}^{{\sf 8c}}_{{\sf
BBQ-K}}$ is symmetric in parameter space the following conclusions will also
hold for the AFM CSL after inversion of the spin direction on one of the two
sublattices on the honeycomb lattice.

We first compare energies in absence of Kitaev interactions at
$\theta/\pi=0.5$. The nearest-neighbor energy levels, exemplified for only
“bright green” spins $\boldsymbol{\sigma}_{\rm g}$ (see Table 1), are given as

$\displaystyle(\boldsymbol{\sigma}_{\rm g},\boldsymbol{\sigma}_{\rm g}):\
E=J_{1}+J_{2}\,,$ (20) $\displaystyle(\boldsymbol{\sigma}_{\rm
g},\bar{\boldsymbol{\sigma}}_{\rm g}):\ E=-J_{1}+J_{2}\,,$ (21)
$\displaystyle(\boldsymbol{\sigma}_{\rm g},\bar{\boldsymbol{\sigma}}_{\rm
r}),\ (\boldsymbol{\sigma}_{\rm g},\bar{\boldsymbol{\sigma}}_{\rm y}),\
(\boldsymbol{\sigma}_{\rm g},\bar{\boldsymbol{\sigma}}_{\rm b}):\
E=\frac{1}{3}J_{1}+\frac{1}{9}J_{2}\,,$ (22)
$\displaystyle(\boldsymbol{\sigma}_{\rm g},\boldsymbol{\sigma}_{\rm r}),\
(\boldsymbol{\sigma}_{\rm g},\boldsymbol{\sigma}_{\rm y}),\
(\boldsymbol{\sigma}_{\rm g},\boldsymbol{\sigma}_{\rm b}):\
E=-\frac{1}{3}J_{1}+\frac{1}{9}J_{2}\,,$ (23)

where the lowest-energy configurations in Eq. (23) exclude nearest neighbors
of same color (e.g., “bright green”–“bright green”) and different contrast
(e.g., “bright green”–“dark red”).

The degeneracy of the lowest-energy states can be lifted by introducing the
bond-dependent Kitaev interactions, giving the ground state energy

$\displaystyle E_{\rm GS}=-\frac{1}{3}J_{1}+\frac{1}{9}J_{2}-\frac{1}{3}K\,.$
(24)

Intriguingly, while $K<0$ gives only one solution for Eq. (24) inducing FM
3$Q$ chiral order, the case of $K>0$ allows for multiple spin pairs
$(\boldsymbol{\sigma}_{\rm g},\boldsymbol{\sigma}_{\rm y})$ and
$(\boldsymbol{\sigma}_{\rm g},\boldsymbol{\sigma}_{\rm r})$ on the $x$ bonds,
$(\boldsymbol{\sigma}_{\rm g},\boldsymbol{\sigma}_{\rm y})$ and
$(\boldsymbol{\sigma}_{\rm g},\boldsymbol{\sigma}_{\rm b})$ on the $y$ bonds,
and $(\boldsymbol{\sigma}_{\rm g},\boldsymbol{\sigma}_{\rm b})$ and
$(\boldsymbol{\sigma}_{\rm g},\boldsymbol{\sigma}_{\rm r})$ on the $z$ bonds.
Since the Kitaev interactions couple spin components on each bond differently,
the allowed color pairs in the ground state will be bond dependent.

Figure 2:  Macroscopic degeneracy in the FM chiral spin liquid (FM CSL). (a)
Constraints on allowed color pairs on bonds for “bright” spins and (b) a
corresponding representative spin configuration in the FM CSL ground-state
manifold. Corresponding figures for “dark” spins are shown in (c) and (d). The
definitions of spin orientations and colors are shown in the insets of (b) and
(d) (see also Table 1).

By doing the same energy analysis for the other “bright” spins
($\boldsymbol{\sigma}_{\rm y}$, $\boldsymbol{\sigma}_{\rm b}$,
$\boldsymbol{\sigma}_{\rm r}$), we obtain the complete set of constraints on
allowed color pairs on bonds:

$\displaystyle(\boldsymbol{\sigma}_{\rm g},\boldsymbol{\sigma}_{\rm
y}),(\boldsymbol{\sigma}_{\rm g},\boldsymbol{\sigma}_{\rm
r}),(\boldsymbol{\sigma}_{\rm b},\boldsymbol{\sigma}_{\rm
y}),(\boldsymbol{\sigma}_{\rm b},\boldsymbol{\sigma}_{\rm r}),$ (25)

on the $x$ bonds,

$\displaystyle(\boldsymbol{\sigma}_{\rm g},\boldsymbol{\sigma}_{\rm
y}),(\boldsymbol{\sigma}_{\rm g},\boldsymbol{\sigma}_{\rm
b}),(\boldsymbol{\sigma}_{\rm r},\boldsymbol{\sigma}_{\rm
y}),(\boldsymbol{\sigma}_{\rm r},\boldsymbol{\sigma}_{\rm b}),$ (26)

on the $y$ bonds, and

$\displaystyle(\boldsymbol{\sigma}_{\rm g},\boldsymbol{\sigma}_{\rm
b}),(\boldsymbol{\sigma}_{\rm g},\boldsymbol{\sigma}_{\rm
r}),(\boldsymbol{\sigma}_{\rm y},\boldsymbol{\sigma}_{\rm
b}),(\boldsymbol{\sigma}_{\rm y},\boldsymbol{\sigma}_{\rm r}),$ (27)

on the $z$ bonds. The same constraints apply for all “dark” colored spins. We
visualize these bond constraints for “bright” and “dark” spins in Figs. 2(a)
and 2(c), respectively.

Every spin in the CSL ground state adheres to these bond constraints, as one
can verify by examining Figs. 2(b) and 2(d), which depict representative spin
configurations in the FM CSL ground-state manifold with "bright" and "dark"
spin sectors, respectively. However, since every spin can choose between two
different colored spins as their neighbors, the bond constraints are not
strong enough to enforce long-range order.

Figure 3:  All allowed arrangements of “bright” spins on a single hexagon in
the FM CSL state, with definitions for colors given in Table 1. Configurations
are arranged in pairs which are related to each other by a six-site hexagon
cluster update following Eq. (28). Definitions of site indices and labels for
the Kitaev bonds are shown on the top left hexagon. Hexagons are shaded in
light cyan, pink, and white, indicating the value of scalar spin chirality
$\kappa_{p}=-8/(3\sqrt{3})$, $16/(3\sqrt{3})$, and $0$, respectively. All
states exhibit $W_{p}=1$ for the semiclassical analog of the
$\mathbb{Z}_{2}$-flux operator. Definitions of $\kappa_{p}$ and $W_{p}$ are
given in Appendix B.

The resulting extensive ground state manifold induces nonzero residual entropy
in the CSL, which can be understood from the fact that two spin configurations
at one hexagon can always be transformed into each other by a local hexagon
cluster update:

${\bf S}^{\rm new}_{p,j}=R^{\alpha}_{j}\ {\bf S}_{p,j}\,,$ (28)

where $R^{\alpha}_{j}$ is the rotation matrix of angle $\pi$ about the axis
$\alpha=x,y,z$; $\alpha$ is the site-dependent label for the Kitaev bond
pointing outwards the hexagon $p$ at site $j=0,\cdots,5$ (see definition on
site labels and $\alpha$ at top left in Fig. 3). In Fig. 3 we show all allowed
spin configurations on a single hexagon in the “bright” manifold of the FM
CSL, which are related by such a six-site hexagon update. Since the CSL allows
always for two spin configurations per hexagon, we obtain the normalized
residual entropy in the ground state

$\frac{S(T\to 0)}{N_{\rm S}}=\frac{1}{N_{\rm S}}\log{2^{N_{\rm
h}}}=\frac{1}{2}\log{2}\,,$ (29)

where $N_{\rm h}=N_{\rm S}/2$ is the number of hexagons on the whole lattice.
This is $1/6$ of the total entropy per site, $\log{8}$.

The set of all possible states on a hexagon in Fig. 3 will also generate a
nonzero scalar spin chirality $\kappa_{p}$, which is defined in Eq. (61) of
Appendix B. The $64$ states per hexagon form three groups with different
values of $\kappa_{p}=-8/(3\sqrt{3})$ (cyan), $16/(3\sqrt{3})$ (pink), and $0$
(white) for $32$, $4$, and $28$ states respectively, which results in the
averaged nonzero value of

$\displaystyle\kappa$ $\displaystyle=\frac{1}{64}\sum_{p}\kappa_{p}$ (30)
$\displaystyle=\frac{1}{64}\left(-\frac{8}{3\sqrt{3}}\times
32+\frac{16}{3\sqrt{3}}\times 4+0\times 28\right)$
$\displaystyle=-\frac{1}{\sqrt{3}}\approx-0.577\,.$

Furthermore, we find that the semiclassical analog of the
$\mathbb{Z}_{2}$-flux operator, $W_{p}$, defined in Eq. (66), gives $W_{p}=1$
for every spin configuration in Fig. 3. We note that $W_{p}=-1$ for the AFM
CSL, offering a way to distinguish the two spin liquids from each other.

Besides the residual entropy and nonzero scalar spin chirality, the bond-
dependent spin constraints in the CSL ground state show extremely short-ranged
spin correlations. As visualized in Fig. 2, nearest-neighbor color pairs are
always of same contrast (“bright”–“bright” or “dark”–“dark”), but never of the
same color. By comparison of spin components, as given in Table 1, the
correlations for all allowed nearest neighbors (distance $\ell=1$) in the
ground-state are therefore always exactly

$D(\ell=1)=\boldsymbol{\sigma}_{i}\cdot\boldsymbol{\sigma}_{j}=-\frac{1}{3}\,.$
(31)

Yet again, by considering the same color constraints, one finds that spin
configurations for next-nearest neighbors and beyond are all equally allowed
and give on average

$D(\ell>1)=\boldsymbol{\sigma}_{i}\cdot\boldsymbol{\sigma}_{k}=3\times\left(-\frac{1}{3}\right)+1\times
1=0\,,$ (32)

where the contribution $1\times 1$ comes from a pair of the same color, which
is forbidden for nearest neighbors, but allowed for further neighbors.
Interestingly, such extremely short-range correlations coincide with
analytical predictions for the $S=1$ Kitaev SL in Ref. [57].

####  III.2.2 Finite-temperature properties of eight-color model

(a) chiral order (b) chiral spin liquid (CSL)

Figure 4:  Comparison of thermodynamic quantities between the FM 3$Q$ chiral
($\phi/\pi=0.2$, $\theta/\pi=0.9$) and AFM 3$Q$ chiral ($\phi/\pi=0.8$,
$\theta/\pi=0.1$) ordered phases (left column), and the FM CSL
($\phi/\pi=0.35$, $\theta/\pi=0.25$) and AFM CSL ($\phi/\pi=0.65$,
$\theta/\pi=0.75$) phases (right column) for the eight-color model given by
$\mathcal{H}^{{\sf 8c}}_{{\sf BBQ-K}}$ in Eq. (15) [see the phase diagram in
Fig. 1(c)]. Shown are the temperature dependencies of (a)–(b) the specific
heat per site, $C/N_{\rm S}$ [Eq. (53)], (c)–(d) the entropy per site,
$S/N_{\rm S}$ [Eq. (54)], and (e)–(f) the absolute value of the scalar spin
chirality, $|\kappa|$ [Eq. (60)]. (g) shows the spin structure factor $S_{\rm
S}({\bf Q})$ [Eqs. (55) and (59)] for the chiral ordered phases at
characteristic ordering vectors indicated in (i) and (k). (h) shows the
semiclassical analog of the $\mathbb{Z}_{2}$-flux operator [Eq. (66)] for the
CSL phases. The momentum resolved spin structure factors, $S_{\rm S}({\bf
q})$, are shown for (i) the FM 3$Q$ chiral ordered phase, (j) the FM CSL
phase, (k) the AFM 3$Q$ chiral ordered phase, and (l) the AFM CSL phase, for
three different temperatures indicated by black arrows in (a) and (b). Data
were obtained from classical MC simulations on clusters of size $L=12,24$, and
$48$ ($N_{\rm S}=288,1152$, and $4608$). $S_{\rm S}({\bf q})$ are shown for
$L=48$, except for (i) and (k) at $T=0.01$, where $L=12$.

In Fig. 4 we show finite-temperature results from classical MC simulations of
$\mathcal{H}^{{\sf 8c}}_{{\sf BBQ-K}}$ in Eq. (15), for chiral ordered phases
(FM 3$Q$ chiral at $\phi/\pi=0.2$ and $\theta/\pi=0.9$, and AFM 3$Q$ chiral at
$\phi/\pi=0.8$ and $\theta/\pi=0.1$), in the left column, and chiral liquid
phases (FM CSL at $\phi/\pi=0.35$ and $\theta/\pi=0.25$, and AFM CSL at
$\phi/\pi=0.65$ and $\theta/\pi=0.75$) in the right column for finite-size
systems of linear dimension $L=12,24$, and $48$. Since the Hamiltonian is
symmetric in parameter space, the FM and AFM states are related to each other
by inversion of the spins on one sublattice of the honeycomb lattice.
Therefore, energy related observables, such as the specific heat $C$ [Eq.
(53)] and the entropy $S$ [Eq. (54)] are identical for the FM and AFM states
(here shown only once), while magnetic observables, like the scalar spin
chirality $\kappa$ [Eq. (60)], the spin structure factor $S_{\rm{S}}({\bf q})$
[Eqs. (55) and (59)], and the semiclassical analog of the $W_{p}$ flux
operator [Eq. (66)] are distinct.

The specific heat, $C$ in the first row of Fig. 4, exhibits a sharp
singularity at (a) $T_{c}=0.186(5)$ and (b) $T_{c}=0.192(5)$ 222 Values in
parentheses represent the statistical error in the last digit, determined by
the standard deviation among five independent simulation runs.. At these
transition temperatures, the system spontaneously breaks its chirality,
selecting a configuration where spins are either all bright or all dark in the
FM 3$Q$ chiral ordered and FM CSL states [see Eq. (23) and discussion in Sec.
III.2.1]. Due to the symmetry of the Hamiltonian in parameter space the AFM
counterparts of these states similarly break chirality. However, in the AFM
versions, there is a mixture of bright and dark spins due to the inversion of
spins on one sublattice of the honeycomb lattice. We analyzed the critical
behavior of these transitions and observed that they do not conform to the
conventional Ising universality class. We leave a more detailed analysis for
future work.

Even though, the spontaneous breaking of chirality seems identical, the left
column of Fig. 4 shows a transition into a magnetically ordered state, while
the right column corresponds to a transition into the disordered CSL state.
This can be directly confirmed in the second row of Fig. 4, which shows the
entropy, $S$ [see Eq. (54)], as obtained from numerical integration of the
specific heat over the available temperature range $10^{-2}\leq T\leq 10$. The
total entropy of the system, $\log{8}$, is released by passing through the
phase transition at $T_{c}$, resulting in zero entropy $S(T\to 0)/N_{\rm S}=0$
in the chiral ordered states as shown in Fig. 4(c). In contrast, the CSLs in
Fig. 4(d) show a residual entropy of

$\frac{S(T\to 0)}{N_{\rm S}}=0.347(1)\approx\frac{1}{2}\log{2}\,,\\\ $ (33)

which matches with the analytically predicted value in Eq. (29), and is direct
numerical evidence of a disordered ground state with extensive degeneracy.

Furthermore, in Figs. 4(e) and 4(f) we show the absolute value of the scalar
spin chirality, $|\kappa|$ [see Eq. (60)], for the ordered and disordered
chiral states. In both cases, $|\kappa|$ scales to zero above $T_{c}$, and
becomes nonzero below $T_{c}$. At low temperatures, we measure for the FM and
AFM 3$Q$ chiral ordered states the following values

$\displaystyle|\kappa|^{\rm FM\ 3{\it Q}\ chiral}$
$\displaystyle=1.5396(1)\approx\frac{8}{3\sqrt{3}}\,,$ (34)
$\displaystyle|\kappa|^{\rm AFM\ 3{\it Q}\ chiral}$
$\displaystyle=3.0792(1)\approx\frac{16}{3\sqrt{3}}\,,$ (35)

and for the FM and AFM CSL states commonly

$|\kappa|^{\rm CSL}=0.5774(1)\approx\frac{1}{\sqrt{3}}\,.$ (36)

Yet again, the result for the CSL matches very well with our analytical
estimate in Eq. (30). Furthermore, as shown in Fig. 4(h), the semiclassical
analog of the $\mathbb{Z}_{2}$-flux operator reaches perfectly $W_{p}=1$ in
the FM CSL, and $W_{p}=-1$ in the AFM CSL, as expected analytically.

Characteristic magnetic signatures, which clearly distinguish ordered from
disordered states, can be found in the dipolar equal-time structure factor
$S_{\rm{S}}({\bf q})$ [see Eqs. (55) and (59)]. We show $S_{\rm{S}}({\bf q})$
at three different temperatures indicated by black arrows in Fig. 4(a), for
the FM and AFM 3$Q$ chiral ordered phases in Figs. 4(i) and 4(k),
respectively. By reducing the temperature, the spectral weight accumulates
around the M points in the Brillouin zone. Below $T_{c}$, Bragg peaks develop
at the M1, M2, and M3 points in the AFM 3$Q$ chiral state. Additionally, Bragg
peaks appear at the M${}^{\prime}_{1}$, M${}^{\prime}_{2}$, and
M${}^{\prime}_{3}$ points in the FM 3$Q$ chiral state. These clearly identify
both states as triple-$q$ ordered states. The intensities of selected Bragg
peaks are explicitly shown as functions of temperature in Fig. 4(g).

Meanwhile, $S_{\rm{S}}({\bf q})$ for the CSL phases behave very differently.
In Figs. 4(j) and 4(l) we show the $S_{\rm{S}}({\bf q})$ for the FM and AFM
CSL states, respectively, at three different temperatures indicated by black
arrows in Fig. 4(b). The signal is very diffuse and does not show any
qualitative difference below and above $T_{c}$. In fact, the scattering
function is reminiscent of the $S_{\rm{S}}({\bf q})$ for the Kitaev SL [89],
showing very diffuse scattering around the Brillouin zone edge in the FM CSL,
and accumulation of intensity around the $\Gamma$ point in the AFM CSL. This
diffuse scattering signature and the absence of any Bragg peaks are the direct
evidence of a disordered state below $T_{c}$.

Figure 5:  Semilogarithmic plot of spin-spin correlations $|D(\ell)|$ [see Eq.
(37)] for the FM CSL case at five different temperatures. Correlations show an
exponential decay above the critical temperature $T_{c}=0.192(5)$ [peak in
Fig. 4(b)]. Below $T_{\rm c}$, correlations show a large value for nearest-
neighbor contributions, $|D(\ell=1)|=0.3333(1)$, and a negligibly small value
beyond nearest neighbors, $|D(\ell>1)|\lesssim 10^{-4}$. Data were obtained
from classical MC simulations of the eight-color model given by
$\mathcal{H}^{{\sf 8c}}_{{\sf BBQ-K}}$ in Eq. (15), at model parameters
$\phi/\pi=0.35$ and $\theta/\pi=0.25$, for a cluster of size $L=48$ ($N_{\rm
S}=4608$). Correlations are shown as function of the Manhatten distance
$\ell$, and were obtained after averaging over symmetrically equivalent paths
along zigzag chains on the honeycomb lattice.

Since momentum-resolved spin-spin correlations of the CSL seem very similar to
those in the Kitaev SL, we investigate correlations also in real space. Figure
5 shows the MC data of the spin-spin correlations

$D(\ell)=\langle\boldsymbol{\sigma}_{0}\cdot\boldsymbol{\sigma}_{\ell}\rangle\,,$
(37)

in the FM CSL at $\phi/\pi=0.35$ and $\theta/\pi=0.25$ measured at five
different temperatures. The correlations are averaged over symmetrically
equivalent paths and are plotted semilogarithmically, where $\ell$ corresponds
to the Manhatten distance [90] along zigzag chains on the honeycomb lattice.
For temperatures above $T_{c}$, correlations decay exponentially, with
increasing correlation length upon reduction of $T$. Below $T_{c}$, the
behavior changes dramatically and exhibits strong nearest-neighbor
correlations with $|D(\ell=1)|=0.3333(1)$ at the lowest measured temperatures
$T=0.01$, while correlations beyond the nearest neighbors become negligibly
small with $|D(\ell>1)|\lesssim 10^{-4}$. These results confirm our analytical
prediction in Eqs. (31) and (32) of extremely short-range correlations in the
CSL.

In short summary of this subsection, our investigation focused on a simplified
model, the eight-color model, $\mathcal{H}^{{\sf 8c}}_{{\sf BBQ-K}}$ in Eq.
(15), with the aim of unraveling the unique nature and nontrivial
characteristics of the eight-color CSL. We discovered that the CSL is
characterized by bond-dependent spin constraints, as visualized in Fig. 2,
giving rise to exotic properties such as residual entropy, nonzero scalar spin
chirality, nonzero $\mathbb{Z}_{2}$-flux order parameter, and extremely short-
range spin correlations. We have successfully verified the existence of all
these properties through unbiased numerical MC simulations.

###  III.3 Analysis of $\mathbb{CP}^{1}$ model

In the previous section, we investigated the eight-color model, as a
simplification of the original BBQ-Kitaev model,
$\mathcal{H}^{\mathcal{S}}_{{\sf BBQ-K}}$ in Eq. (1), by discretising the
continuous spin space of an $S=1$ moment to only eight allowed, noncoplanar
states. This helped us to transparently understand the nature and origin of
the exotic CSL phase in a simple analytical way.

In this section, we are going to partially relax this constraint, and allow
dipolar spin components to take continuous values, namely being a classical
Heisenberg vector described in the order parameter space $\mathbb{CP}^{1}$.
Due to the enlarged degree of freedom in spin space, the model offers a
variety of additional and nontrivial phases. Understanding these phases will
help us in bridging our comprehension from the simplified eight-color model to
the more complex $\mathbb{CP}^{2}$ model.

####  III.3.1 Ground-state properties of $\mathbb{CP}^{1}$ model

Figure 6:  Competing phases of $\mathcal{H}^{\mathcal{S}}_{{\sf BBQ-K}}$ in
Eq. (1) for classical Heisenberg spins ($\mathbb{CP}^{1}$ model), obtained by
variational energy minimization for a finite-size cluster of linear dimension
$L=120$ ($N_{\rm S}=28\ 800$ spins) under periodic boundary conditions. (a)
Zoom-in of the ground-state phase diagram in Fig. 1(b) to the region for
$0.0\leq\phi/\pi\leq 1.0$. (b) Normalized energy and its derivative, $E/N_{\rm
S}$ and $\partial(E/N_{\rm S})/\partial\phi$, (top), and the structure factor
for dipole moments, $S_{\rm S}({\bf q})$ [Eqs. (55) and (59)], at high-
symmetry momenta (bottom) at $\theta/\pi=0.06$ [red line in (a)]. The solid
and dashed lines represent first- and second-order phase transitions,
respectively. (c) Energy difference, $\Delta E/N_{\rm S}$, measured from the
ground state to the metastable FM CSL and 1D-1 phases at $\theta/\pi=0.06$.
The lower insets show spin configurations in each ordered phase, while the
upper ones represent metastable states of the FM CSL.

In Fig. 6(a) we show a zoom-in of the ground-state phase diagram in Fig. 1(b)
for positive $J_{2}$ $(0.0<\phi/\pi<1.0)$, where the model is mostly
frustrated. Between FM/AFM ordered and FM/AFM 3$Q$ chiral phases the model
stabilizes NC phases, and canted planar (CPL) phases over a wide range of
model parameters. The canting angle of spins in those phases depends on the
model parameters. Notably, specific choices of $\phi$ allow us to stabilize
disordered CSLs and dimensionally-reduced 1D phases in the ground state, as
indicated by the black dotted lines in Fig. 6(a) 333 We observed several
unconventional phases near the two SU(3) points at
$(\theta/\pi,\phi/\pi)=(0.5,0.25)$ and $(0.5,0.5)$, which we reserve for a
more in-depth discussion in the future..

To see these phase competitions in more detail, in Fig. 6(b) we concentrate on
the lower part of the phase diagram at $\theta/\pi=0.06$ [red line in Fig.
6(a)], and explicitly show the normalized energy $E/N_{\rm S}$ and its
derivative $\partial(E/N_{\rm S})/\partial\phi$ (top), together with the
$S_{\rm S}({\bf q})$ [Eqs. (55) and (59)] at relevant ordering vectors K′, M1,
M2, M3, and M${}^{\prime}_{3}$ (bottom). The transition from CPL-1 to AFM 3$Q$
chiral phase at $\phi/\pi\approx 0.611$ is of first order (solid black line),
while transitions from AFM to NC-1 at $\phi/\pi\approx 0.148$ and NC-1 to
CPL-1 at $\phi/\pi=0.5$ are of second order (dashed black lines), as seen by a
jump and kinks in $\partial(E/N_{\rm S})/\partial\phi$, respectively. The
Bragg peaks at the K′ points monotonically reduce their intensity in the NC-1
phase as $\phi$ increases, while additional ordering vectors at the M2 and
M${}^{\prime}_{3}$ points increase. In the CPL-1 phase, intensities at the M3
and M1 points become also nonzero, while the AFM 3$Q$ chiral state shows an
equal weight at the M1, M2, and M3 points, forming the triple-$q$ order, as
also obtained in the eight-color model in Fig. 4(k).

In Fig. 6(c) we show the energy difference $\Delta E/N_{\rm S}$ between the FM
CSL (yellow circles) and the 1D-1 state (blue circles) relative to the ground
state for $\theta/\pi=0.06$. The lower insets illustrate the spin
configurations within each ordered phase, while the upper insets represent
metastable states of the FM CSL. Upon increasing $\phi/\pi$ above $\approx
0.148$, spins in the ground state start to cant away from the collinear AFM
arrangement to establish the NC-1 order, with a four-site magnetic unit cell
[see Fig. 1(d)]. The canting angle of spins in the NC-1 phase is $\phi$
dependent, and the spontaneous selection of a particular direction for this
canting angle induces a single-$q$ order with Bragg peaks at the M2 points in
the first Brillouin zone, supplemented by other peaks in the extended
Brillouin zone [see Fig. 6(b)].

Upon tuning the model to $\phi/\pi=(1/\pi)\arctan(3/2)$ [see Eq. (16)], the
spins achieve isotropic alignment across the lattice, by pointing precisely
towards the corners of a unit cube. At this critical point, the system allows
an extensive number of spin configurations to form the CSL ground state,
characterized by all properties as discussed on the eight-color model in Sec.
III.2.1. In the vicinity of this critical point, we find a CSL like state as a
metastable state. Spins in this metastable CSL are slightly canted away from
the (111) directions as shown in the insets of Fig. 6(c), while still
retaining the underlying nature of the CSL with macroscopic degeneracy.

Upon further increasing $\phi$, spins in the NC-1 ground state continue to
cant until they achieve complete coplanarity, precisely at the boundary
between the NC-1 and CPL-1 phases, at $\phi/\pi=0.5$. At this specific point,
nearest-neighbor spins are mutually orthogonal to each other, and effectively
cancel nearest-neighbor bilinear and biquadratic interactions. Coplanar spins
globally select either the $xy$, $yz$, or $zx$ plane [the inset of Fig. 6(c)
shows the case of the $xy$ plane], by setting one of their spin components to
zero. Therefore the dominating Kitaev interactions for one particular bond-
type will be zero, and effectively decouple individual zigzag chains
throughout the whole lattice. The system effectively forms a dimensionally-
reduced 1D phase. Increasing $\phi$ beyond this point will again impose
canting of spins, but this time with an enlarged eight-site magnetic unit cell
forming the magnetic order of the CPL-1 state [see Fig. 1(d)].

####  III.3.2 Finite-temperature properties of $\mathbb{CP}^{1}$ model

Figure 7:  Finite-temperature properties of the FM CSL phase for the
$\mathbb{CP}^{1}$ model given by $\mathcal{H}^{\mathcal{S}}_{{\sf BBQ-K}}$ in
Eq. (1) at $\phi/\pi=(1/\pi)\arctan(3/2)$ and $\theta/\pi=0.06$. Shown are
temperature dependencies of (a) the specific heat per site, $C/N_{\rm S}$ [Eq.
(53)], (b) the absolute value of the scalar spin chirality $|\kappa|$ [Eq.
(60)], and (c) the semiclassical analog of the $\mathbb{Z}_{2}$-flux operator,
$W_{p}$ [Eq. (66)]. MC simulations were performed for $L=6,12,24$, and $48$
($N_{\rm S}=72,288,1152$, and $4608$). (d) shows the spin structure factor
$S_{\rm S}({\bf q})$ [Eqs. (55) and (59)] at temperatures indicated by black
arrows in (a) for $L=24$.

Let us now turn to the thermodynamic properties of the CSL in the
$\mathbb{CP}^{1}$ model. In Fig. 7, we show results from classical MC
simulations of $\mathcal{H}^{\mathcal{S}}_{{\sf BBQ-K}}$ in Eq (1), at model
parameters $\theta/\pi=0.06$ and $\phi/\pi=(1/\pi)\arctan(3/2)$, which
stabilize the FM CSL in the ground state [see Eq. (16)].

In Figs. 7(a)–7(c) we respectively show the specific heat $C/N_{\rm S}$ [Eq.
(53)], the scalar spin chirality $|\kappa|$ [Eq. (60)], and the semiclassical
analog of the $\mathbb{Z}_{2}$-flux operator, $W_{p}$ [Eq. (66)], over a wide
range of temperatures. $C/N_{\rm S}$ shows a clear singularity at
$T_{c}=0.0193(6)$, which scales weakly with system size, and has a shape which
looks somewhat similar to the observed singularity in the eight-color model in
Fig. 4(b).

For $T\to 0$, we observe $C/N_{\rm S}\to 1$, as expected for classical
Heisenberg spins in the absence of soft mode excitations. We note that the
acceptance ratio for single-spin flip updates at such low temperatures becomes
strongly suppressed (see Appendix D), suggesting that continuous local spin
motion is not able to adiabatically connect different ground-state
configurations. However, where the single-spin flip fails, the hexagon update
becomes successful, demonstrating that in the CSL phase, decorrelation
primarily occurs through cluster updates [see Eq. (28)].

The scalar spin chirality $|\kappa|$ scales to zero in the paramagnetic
regime, and takes nonzero values below $T_{c}$, as observed in the eight-color
model in Fig. 4(f). However, in comparison, $|\kappa|$ does not rapidly
saturate, but rather monotonically reaches a value slightly smaller than in
the eight-color model [see Eq. (36)]. We note that the asymptotic value of
$|\kappa|$ becomes $\theta$ dependent (not shown here). However, variational
energy minimization confirms that at $T=0$, spins precisely align along the
(111) direction, as also illustrated in the lower insets of Fig. 6(c). This
seemingly contradictory observation implies that different configurations of
the CSL ground state manifold in Fig. 3 are selected with varying weights in
the $\mathbb{CP}^{1}$ model, resulting in a reduction of $\kappa$ at finite
$T$. At the lowest temperatures, the $\mathbb{Z}_{2}$-flux operator takes a
value $W_{p}\approx 1$, as observed in the eight-color model.

In Fig. 7(d) we show the spin structure factor $S_{\rm S}({\bf q})$ at three
different temperatures indicated by black arrows in Fig. 7(a). The signal is
diffuse across the entire range of temperatures and closely resembles the spin
structure factor of the eight-color FM CSL [see Fig. 4(j)] with an enlarged
concentration of intensity around the K’ points and, almost invisible, some
additional weak intensity accumulation around the M points in the first
Brillouin zone at $T=0.022$.

Thermodynamic observables for the CSL in both the eight-color model and the
$\mathbb{CP}^{1}$ model match remarkably well. This alignment provides strong
evidence that the eight-color CSL, with all of its physical properties
analyzed in Sec. III.2.1, remains preserved even after relaxing the spin
degree of freedom to $\mathbb{CP}^{1}$. Given that the CSL persists as a
metastable state in the vicinity of its optimal model parameters given by Eq.
(16) [see Fig. 6(c)], we anticipate its survival for an enlarged region of
model parameters at finite temperatures due to the substantial entropy
associated with its macroscopic degeneracy.

###  III.4 Analysis of $\mathbb{CP}^{2}$ model

The original BBQ-Kitaev model contains both dipole and quadrupole fluctuations
of $S=1$ magnetic moments, which are correctly described in the spin space
$\mathbb{CP}^{2}$. In the last two subsections we simplified the problem by
restricting the spin degree of freedom to extract the essential properties of
the CSL state. In the present subsection we come back to the original model by
reintroducing all allowed degrees of freedom for $S=1$ moments. We will not
only see that the CSL survives at finite temperatures even in the presence of
quadrupolar degree of freedom, but also quantitatively examine the impact of
local quantum fluctuations on the ground-state phase diagram. Some results in
this subsection have already been presented in our earlier work in Ref. [74],
but we include some of them here to make the discussion and comparison to
$\mathbb{CP}^{1}$ and eight-color models self-contained.

####  III.4.1 Ground-state properties of $\mathbb{CP}^{2}$ model

Figure 8:  Competing phases of $\mathcal{H}^{\mathcal{A}}_{{\sf BBQ-K}}$ in
Eq. (7) for $S=1$ magnetic moments ($\mathbb{CP}^{2}$ model), obtained by
variational energy minimization for a finite-size cluster of linear dimension
$L=120$ ($N_{\rm S}=28\ 800$ spins) under periodic boundary conditions. (a)
Averaged spin norm $|{\bf S}|$ [see Eq. (10)] in the region
$0.2\leq\phi/\pi\leq 0.5$. The solid and dashed lines represent first- and
second-order phase transitions, respectively. (b) The normalized energy
$E/N_{\rm S}$ and its derivative $d(E/N_{\rm S})/d\phi$ (top panel), norms for
dipole moment, $|{\bf S}|$, and quadrupole moment, $|{\bf Q}|$ [see Eq. (11)]
(middle panel), and the structure factors for dipole moments, $S_{\rm{S}}({\bf
q})$ [Eqs. (55) and (56)], and quadrupole moments, $S_{\rm{Q}}({\bf q})$ [Eqs.
(55) and (57)] (bottom panel), along $\theta/\pi=0.06$ [red line in (a)]. (c)
Energy difference $\Delta E/N_{\rm S}$, measured from the NC-1 state. The
lower insets show spin configurations in each ordered phase, while the upper
ones resemble metastable states of the FM CSL.

In Fig. 8(a), we show the spin-dipole norm $|{\bf S}|$ in the frustrated
region of the phase diagram in Fig. 1(a). In most of the phases, the spin norm
is $|{\bf S}|\neq 1$, indicating the presence of contribution from onsite
quadrupole moments. The onsite dipole and quadrupole moments mix in nontrivial
ways to alleviate frustration, leading to quantitative changes in phases and
their boundaries, in contrast to the $\mathbb{CP}^{1}$ model shown in Fig.
6(a).

The purely quadrupolar SO phase at $\theta/\pi=0.5$ [see Fig. 1(d)],
immediately changes into a q1D ordered state by introducing the Kitaev
interactions. For both $\theta/\pi>0.5$ and $\theta/\pi<0.5$, the length of
dipolar spin moments monotonically increases and behaves essentially the same
as described for the 1D phase in the $\mathbb{CP}^{1}$ model in Sec. III.3.1,
where spins form dominant correlations along isolated 1D zigzag chains. The
key distinction from the $\mathbb{CP}^{1}$ model lies in the influence of
quadrupole correlations, which introduce weak interactions between the chains
at very low temperatures. This leads to the formation of a 2D ordered state
exhibiting one-dimenaional character over a wide region at low and
intermediate temperatures. This characteristic is the reason why we refer to
the state as a quasi-1D state. Significantly, these quadrupole correlations
stabilize the region of the q1D phases within the phase diagram, a stark
contrast to the $\mathbb{CP}^{1}$ model, where the 1D phases remains stable
only along a singular line. However, to this expense, the NC phases are
notably suppressed, along with the neighboring TC phases.

In Fig. 8(b), we focus along the line at $\theta/\pi=0.06$ [red line in Fig.
8(a)], and plot $E/N_{\rm S}$ and $\partial(E/N_{\rm S})/\partial\phi$ (top
panel), the spin dipole norm $|{\bf S}|$ [see Eq. (10)] and quadrupole norm
$|{\bf Q}|$ [see Eq. (11)] (middle panel), and the structure factors for
dipoles, $S_{\rm S}({\bf q})$ [Eqs. (55) and (56)], and for quadrupoles,
$S_{\rm Q}({\bf q})$ [Eqs. (55) and (57)], at relevant ordering vectors
(bottom panel). We identify first-order phase transitions (solid lines) from
NC-1 to TC-1 at $\phi/\pi=0.3085(5)$, and from q1D-1 to AFM 3$Q$ chiral at
$\phi/\pi=0.4015(1)$. Additionally, we find second-order transitions (dashed
lines) from AFM to NC-1 at $\phi/\pi=0.250(1)$, and from TC-1 to q1D-1 at
$\phi/\pi=0.3305(5)$. Notably, our measurements reveal a reduction in the spin
dipole norm $|{\bf S}|<1$ within the NC-1, TC-1, and q1D-1 phases. It is worth
noting that the reduction of $|{\bf S}|$ is $\theta$ dependent; however, in
this parameter region near the Kitaev limit, it is relatively small, with less
than $1$ %. This highlights the significance of dipole moments and justifies
our comparison to the simplified $\mathbb{CP}^{1}$ model for dominant Kitaev
interactions, as done in Sec. III.3.

In Fig. 8(c), we present the energy difference between the metastable FM CSL
and TC-1 states to the NC-1 ordered state. The presence of local quadrupole
moments induces a small, yet finite energy difference to the FM CSL bounded by

$\frac{\Delta E}{N_{\rm S}}=0.000393(4)\,,$ (38)

for $\phi/\pi=0.305(1)$. Consequently, the CSL becomes an excited state rather
than the ground state, with spin configurations no longer perfectly aligned
with the corners of the unit cube, but still retaining the essential
properties of the CSL as a metastable state. In fact, we find that the FM CSL
is a metastable state with small energy gap in a wide region of the phase
diagram, and smoothly merges with the AFM Kitaev SL for $\theta\to 0$ (see the
Supplemental Material in Ref. [74]). The evolution of spin canting angles,
illustrated in the bottom insets, closely resembles the general behavior
discussed in the $\mathbb{CP}^{1}$ model [see Fig. 6(c)]. We note, however, a
noticeable difference: At $\phi/\pi=0.3085(5)$, the TC-1 state [see spin
configurations in lower inset of Fig. 8(c)] becomes the ground state, a
scenario which seems absent in the $\mathbb{CP}^{1}$ model. Furthermore, the
q1D-1 state maintains its status as the ground state for $\phi/\pi>0.3305(5)$,
in contrast to the $\mathbb{CP}^{1}$ model where it remained the ground state
only along a singular line.

####  III.4.2 Finite-temperature properties of $\mathbb{CP}^{2}$ model

Figure 9:  MC simulations of $\mathcal{H}^{\mathcal{A}}_{{\sf BBQ-K}}$ in Eq.
(7), for the $\mathbb{CP}^{2}$ model at $\phi/\pi=0.305$ and $\theta/\pi=0.06$
reveal the presence of the FM CSL at finite temperatures. Shown are
temperature dependencies of (a) the specific heat per site, $C/N_{\rm S}$ [Eq.
(53)], (b) the absolute value of the scaler spin chirality $|\kappa|$ [Eq.
(60)], and (c) the semiclassical analog of the $\mathbb{Z}_{2}$-flux operator,
$W_{p}$ [Eq. (66)]. Simulations were done for $L=6,12,24$, and $48$ ($N_{\rm
S}=72,288,1152,$ and $4608$). (d) and (e) show the spin structure factor for
dipoles, $S_{\rm S}({\bf q})$, and quadrupoles $S_{\rm Q}({\bf q})$ [Eqs.
(55)–(57)], respectively, at temperatures indicated by black arrows in (a) for
$L=24$ at $T=0.01$ and $T=0.003$, while for $L=12$ at $T=0.001$. The green
path in the left panel of (d) corresponds to the path used for the dynamical
structure factor in Fig 10.

The extensive degeneracy in the FM CSL and the small but finite energy gap to
the ground state suggest that the system undergoes an entropy-driven phase
transition from NC-1 order to the FM CSL at finite temperature. To confirm
this, we show finite-temperature MC simulations for
$\mathcal{H}^{\mathcal{A}}_{{\sf BBQ-K}}$ [Eq. (7)] in Fig. 9, for model
parameters where the FM CSL shows the smallest energy difference to the ground
state ($\theta/\pi=0.06$ and $\phi/\pi=0.305$). Figures 9(a)–9(d) are
basically the same as those in Ref. [74], but we additionally show the static
structure factors for quadrupoles $S_{\rm Q}({\bf q})$ in Fig. 9(e).

The specific heat $C/N_{\rm S}$ [Eq. (53)] in Fig. 9(a) shows two
singularities at $T_{c1}=0.0055(2)$ and $T_{c2}=0.0017(2)$, which scale weakly
with system size. For $T\to 0$, we observe $C/N_{\rm S}\to 2$, as expected for
$S=1$ magnetic moments described within the spin space of $\mathbb{CP}^{2}$ in
the absence of soft-mode fluctuations [47]. Below $T_{c2}$, characteristic
Bragg peaks emerge in the structure factors at the momentum points,
corresponding to the NC-1 order. Therefore, we associate $T_{c2}$ with a
symmetry-breaking transition into the NC-1 ordered phase [74].

The intermediate phase between $T_{c1}$ and $T_{c2}$ shows all properties of
the FM CSL discussed in the previous subsections, namely a nonzero, size-
independent scalar spin chirality $|\kappa|$ [Eq. (60)], a value of $W_{p}$
[Eq. (66)] which is almost $+1$, and a very diffuse structure factor for spin
dipoles, $S_{\rm S}({\bf q})$. The values of both quantities, $|\kappa|$ and
$W_{p}$, are somewhat reduced compared to the ideal eight-color CSL discussed
in Sec. III.2. This reduction stems from the decrease in spin lengths $|{\bf
S}|$ caused by the presence of small, but finite, quadrupole moments [see Fig.
8(b)] and temperature fluctuations.

The diffuse $S_{\rm S}({\bf q})$ is comparable to the structure factors found
in the eight-color and $\mathbb{CP}^{1}$ models [see Figs. 4(j) and 7(d),
respectively], albeit some additional weak intensity features around the M
points in the first Brillouin zone. The quadrupole structure factor, $S_{\rm
Q}({\bf q})$, also suggests the absence of quadrupolar order, with scattering
patterns closely resembling those for dipoles. This similarity arises because
spin moments are primarily dipolar in nature, and dipole characteristics are
implicitly reflected in the quadrupole components. Taken these observations,
we associate $T_{c1}$ with a discrete chiral symmetry breaking into the FM
CSL, as observed in the eight-color model and $\mathbb{CP}^{1}$ model cases
[see Figs. 4(b) and 7(a), respectively]. It is worth noting that the
intermediate-temperature CSL persists over an extended region around
$\phi/\pi\approx 0.3$ and $\theta/\pi\approx 0.2$, as shown in the
Supplemental Material of Ref. [74].

The scattering in $S_{\rm S}({\bf q})$ at $T=0.01$, above $T_{c1}$, exhibits
characteristics reminiscent of the semiclassical AFM Kitaev SL, as our chosen
model parameters contain strong Kitaev interactions. The intensity exhibits
diffuse accumulation around the Brillouin zone edge, similar to the results
for the “pure” Kitaev model in Appendix E.1, however, with a stronger
accumulation of diffuse intensity around the K’ points. Meanwhile, the
structure factor $S_{\rm Q}({\bf q})$ shows a more uniform intensity
distribution with dominant scattering intensity around the K and K’ points.

All physical observables characterizing the CSL in the $\mathbb{CP}^{2}$ model
closely resemble those measured in the eight-color model (Sec. III.2) and the
$\mathbb{CP}^{1}$ model (Sec. III.3). This remarkable consistency strongly
suggests that the dominant properties of the CSL, as discussed in Sec.
III.2.1, remains robust even after considering all allowed local degrees of
freedom for an $S=1$ moment in the spin space $\mathbb{CP}^{2}$. However, it
is important to emphasize that the CSL no longer represents the ground state
but emerges as an entropically driven state at finite temperature.

###  III.5 Dynamical properties of the chiral spin liquid

In this section we analyze dynamical properties of the CSL in the
$\mathbb{CP}^{2}$ model. We show that the excitations corresponding to the
CSL, in fact, are gapped, consistent with our observation of extreme short-
range correlations in Fig. 5.

####  III.5.1 Excitation spectrum

(a)  (b)  (c)

(d)  (e)  (f)

Figure 10:  Dynamical structure factors [see Eqs. (67)–(69)] for spin-dipole
moments, ${\tilde{S}}_{\rm{S}}({\bf q},\omega)$ (top panels), and spin-
quadrupole moments, $\tilde{S}_{\rm{Q}}({\bf q},\omega)$ (bottom panels) of
$\mathcal{H}^{\mathcal{A}}_{{\sf BBQ-K}}$ in Eq. (7), for the
$\mathbb{CP}^{2}$ model, at $\phi/\pi=0.305$ and $\theta/\pi=0.06$. Results
are shown for (a) and (d) the NC-1 ordered phase at $T=0.001$, (b) and (e) the
FM CSL at $T=0.003$, and (c) and (f) the semiclassical analog of the AFM
Kitaev SL at $T=0.01$. Temperatures are indicated by black arrows in Fig.
9(a). Data were obtained from molecular dynamics simulations (details in
Appendix A) for finite-size clusters of $L=48$ $(N=4608)$. Results are plotted
along the green path in momentum space, as indicated in the left panel of Fig.
9(d).

In Fig. 10, we show the dynamical structure factors for spin-dipole moments,
$\tilde{S}_{\rm{S}}({\bf q},\omega)$ (top panels), and spin-quadrupole
moments, $\tilde{S}_{\rm{Q}}({\bf q},\omega)$ (bottom panels), which are
obtained from molecular dynamics simulations (see details in Appendix A) at
the same model parameters as chosen in Fig. 9. We show the spectrum along the
path in momentum space, as indicated by green lines in the left panel of Fig.
9(d) and compare dynamics between the NC-1 ordered phase at $T=0.001$, the FM
CSL at $T=0.003$, and the semiclassical analog of the AFM Kitaev SL at
$T=0.01$. Their corresponding equal-time structure factors are plotted in
Figs. 9(d) and 9(e).

Figures 10(a) and 10(d) show the dynamical structure factors for the NC-1
ordered phase. For model parameters used here, the energy cost of a single-
spin flip is primarily controlled by the Kitaev bond energies, limiting the
bandwidth of excitations to $\omega/J\approx 2$. Given the magnetic order in
this phase, well-defined spin waves and quadrupole waves emerge in the dipole
and quadrupole channels, respectively. The signal reveals a nontrivial
dispersion, which can be separated into three energy regimes. In the low-
energy regime for $\omega/J\lesssim 0.5$, linearly dispersing Goldstone modes
exist at the magnetic ordering vectors $\Gamma$, M2, and K’ (M’3 is not shown
here), in accordance to the Bragg-peaks seen in the left column of Figs. 9(d)
and 9(e). Moving to higher energy between $0.6\lesssim\omega/J\lesssim 1.9$,
two well-defined gapped bands exist, equally present in spin-dipole as well as
spin-quadrupole channels. At the energy maximum, the spectrum features an
almost flat band, corresponding to quadrupole excitations captured by strong
intensity in $\tilde{S}_{\rm{Q}}({\bf q},\omega)$. The band is almost flat, as
the biquadratic interactions $J_{2}$ are small compared to the dominant Kitaev
interactions $K$. A weak-intensity “shadow” of this excitation is visible in
the dipole channel, coming from mixing between dipole and quadrupole
excitations due to spin-anisotropic interactions.

When the system turns into the FM CSL by increasing the temperature above
$T_{c2}$, noticeable changes occur in the dynamical signatures, as shown in
Figs. 10(b) and 10(e). The once well-defined spin waves at higher energy in
the NC-1 ordered phase form a broad continuum of diffuse excitations with a
gap of approximately $1.2J$ and a band maximum of around $1.9J$. We associate
this continuum with the gapped excitations of the FM CSL. Considering the
simplified eight-color model in Sec. III.2, an elementary excitation out of
the CSL ground-state manifold is comprised of a single-spin flip, violating
the ground-state constraints in Fig. 2. However, such a single-spin flip will
violate the ground-state constraints not only on one but on two bonds, with an
energy cost on each bond of $\frac{2}{3}K$ [see Eq. (24)]. By using
$K=\cos{(0.06\pi)}$, for the model parameters used in this calculation (see
Fig. 9), we estimate a minimum excitation energy of

$\Delta\omega^{\rm CSL}_{\rm min}=2\times\frac{2}{3}K\approx 1.31J\,.$ (39)

This estimate aligns well with the energy gap observed in Fig. 10(b). Above
this energy region, quadrupole excitations do not form a continuum but instead
exhibit well-defined quadrupole waves, with a remnant of one of the two
branches seen in the NC-1 phase in Fig. 10(d). This is rather intriguing,
since the energy-integrated correlations in Figs. 9(d) and 9(e) show a diffuse
structure which is very similar for both $S_{\rm S}({\bf q})$ and $S_{\rm
Q}({\bf q})$. This suggests that excitations in the CSL continuum induce
coherent quadrupole dynamics at finite frequencies.

Furthermore, we observe that the well-defined Goldstone modes, corresponding
to the NC-1 order change into a strong zero-energy flat band and a diffuse
continuum of excitations with a small energy gap of $\sim 0.1J$ and bandwidth
of $\sim 0.5J$. As the zero-energy mode is challenging to distinguish from the
diffuse low-energy continuum, we provide the dynamics for the FM CSL in the
$\mathbb{CP}^{1}$ model in Appendix E.2. Here, dynamics exhibit more
pronounced and sharper bands, allowing us to validate their genuine physical
presence and ensuring that they are not artifacts of the simulations.

In Figs. 10(c) and 10(f), we show results above $T_{c1}$ in the semiclassical
analog of the AFM Kitaev SL. It appears that the three energy regimes observed
in the FM CSL phase merge to form a spin structure factor that spans the full
energy range up to $\omega/J\approx 2$, with strong intensity below $\sim
0.5J$. The quadrupole channel is dominated by an almost flat band at
$\omega/J\approx 2$, which corresponds to the same excitations of quadrupolar
origin as discussed in Figs. 10(d) and 10(e). The observed dynamical
correlations in this finite-temperature phase strongly resemble the dynamics
of the semiclassical Kitaev SL, as shown in the comparison in Appendix E.1.
This similarity supports our interpretation of a phase transition from the
semiclassical Kitaev SL to the CSL state, since chosen model parameters ensure
that Kitaev interactions strongly favor Kitaev SL physics above $T_{c1}$.

####  III.5.2 Energy cross sections

Figure 11:  Energy cross sections of the dynamical structure factors [see Eqs.
(67)–(69)] for the FM CSL phase in the $\mathbb{CP}^{2}$ model of
$\mathcal{H}^{\mathcal{A}}_{{\sf BBQ-K}}$ [see Eq. (7)], as shown in Figs.
10(b) and 10(e). Results are shown in (a)–(c) for dipoles,
$\tilde{S}_{\rm{S}}({\bf q},\omega)$, while in (d)–(f) for quadrupoles,
$\tilde{S}_{\rm{Q}}({\bf q},\omega)$.

In the following, we examine the dynamical signatures in the FM CSL, shown in
Figs. 10(b) and 10(e), at $T=0.003$, in more detail by analyzing energy cross
sections. In Fig. 11 we show dynamical structure factors [see Eqs. (67)–(69)]
for spin-dipole moments, $\tilde{S}_{\rm{S}}({\bf q},\omega)$ (top panels) and
spin-quadrupole moments, $\tilde{S}_{\rm{Q}}({\bf q},\omega)$ (bottom panels)
at particular energy values. The cross section at low energy, $\omega/J=0.02$,
in Figs. 11(a) and 11(d) shows a characteristic high-intensity, diffuse
scattering signal with dominant weight around the K’ points and some
accumulation of weak spectral weight around the M points. These signals are
almost completely captured by the equal-time structure factors $S_{\rm S}({\bf
q})$ and $S_{\rm Q}({\bf q})$, shown in Figs. 9(d) and 9(e) (middle column).

At small, but finite energies, the spectrum shows a diffuse signal extending
up to $\omega/J\approx 0.5$. An energy cross section within this energy range
at $\omega/J=0.26$ reveals complicated scattering features with finite
intensity across the entire Brillouin zone. We observe dominant weight around
the edge of the extended Brillouin zone for dipoles in Fig. 11(b), and
predominant scattering within the first Brillouin zone for quadrupoles in Fig.
11(e).

The energy cross section within the diffuse and almost featureless continuum
for spin dipoles at $\omega/J=1.4$ reveals an almost uniform distribution of
intensity over the entire Brillouin zone in Fig. 11(c). Within this energy
range, where dipoles exhibit a continuum, quadrupole excitations display well-
defined energy bands, as visible in Fig. 10(e). Shown in Fig. 11(f), these
bands manifest as rings around the K’ points in momentum space.

In Appendix E.2, we explicitly compare the dynamical properties of the CSL in
the $\mathbb{CP}^{2}$ model with the results obtained for the
$\mathbb{CP}^{1}$ model. Our findings reveal a qualitative agreement between
the two models, providing additional support for our interpretation that the
spin liquid physics of the CSL is primarily governed by spin dipoles.

##  IV Summary and perspectives

In this paper, we explored the $S=1$ Kitaev model with bilinear-biquadratic
(BBQ) interactions on the honeycomb lattice for SU(3) spin-coherent states in
the spin space $\mathbb{CP}^{2}$, a model which hosts a diversity of exotic
phases, including multiple-$q$ states with nonzero scalar spin chirality,
quasi-one-dimensional coplanar phases, twisted conical phases, and noncoplanar
ordered states. Building upon the theory presented in our prior work, Ref.
[74], our goal was to provide a pedagogical explanation of the nature and
unique properties of a finite-temperature chiral spin liquid (CSL), apparent
in this model. To achieve this, we utilized variational energy minimization
and classical Monte Carlo (MC) with molecular dynamics simulations, using the
recently developed U(3) formalism specifically designed for simulating $S=1$
magnets [47].

By progressively restricting the spin degree of freedom from $\mathbb{CP}^{2}$
to $\mathbb{CP}^{1}$ and eventually to a discrete eight-color model, we
unveiled the significant impact of the local Hilbert space dimension on the
ground state and its excitation properties. The eight-color model, designed to
capture the most intriguing aspects of the CSL, offered an intuitive and
analytical understanding of its physical properties. This model not only
uncovered the origin of the extensive degeneracy in the CSL, but also
clarified the reasons behind the nonzero scalar spin chirality, extreme short-
ranged correlations, and the presence of a Kitaev spin liquid (SL) feature of
$\mathbb{Z}_{2}$-flux order. We derived bond-dependent spin constraints that
unequivocally determine the nature of this eight-color CSL and validated our
findings through unbiased classical MC simulations. Our results convincingly
demonstrated that the dominant properties of the eight-color CSL persist even
in the more complex $\mathbb{CP}^{1}$ and $\mathbb{CP}^{2}$ model cases.
Importantly, the enlarged spin degree of freedom does not alter the
fundamental properties of the CSL but rather restricts its stability in the
ground state. Consequently, the CSL survives as an entropy-driven, robust
phase at finite temperatures in the original $\mathbb{CP}^{2}$ model. This
stability in the more complex $\mathbb{CP}^{2}$ model enabled us to simulate
the dynamical structure factor, which revealed a nontrivial excitation
spectrum characterized by a high-intensity zero-energy mode with a diffuse
continuum of excitations in the low-frequency region. At higher frequencies, a
broad continuum of excitations was observed, which we attributed to the
excitations of the discrete eight-color CSL.

CSLs represent a well-established branch of SLs, wherein the SL state persists
even after the spontaneous breaking of time-reversal symmetry [92, 93]. While
examples of CSLs with exactly known quantum-ground states exist [94, 95],
studies across various models have indicated that the spontaneous breaking of
discrete time-reversal symmetry typically occurs at finite temperatures [96,
97, 98, 99, 100]. The CSL investigated in our present study is a classical SL
that also exhibits the spontaneous breaking of time-reversal symmetry at
finite temperatures, while showing the absence of magnetic order at lower
temperatures. Remarkably, the correlations above and below this phase
transition are strikingly similar [see Figs. 4(i) and 4(l)], and suggest a
potential deeper connection to correlations seen in Kitaev liquids [89].
Despite the critical behavior not conforming to the conventional Ising
universality class, a detailed scaling analysis of this transition remains a
topic for future investigation. It would be intriguing to explore whether the
phase transition in the $\mathbb{CP}^{1}$ and $\mathbb{CP}^{2}$ models
exhibits a similar scaling behavior as observed in the eight-color model, and
whether such a behaviour could also be observed in other lattice models which
host classical CSLs [101]. Furthermore, the internal structure and properties
of the CSL, including the bond-dependent spin constraints illustrated in Fig.
2, suggest the potential existence of a suitable lattice gauge theory [102,
103].

We strongly advocate for simulating systems with higher-order interactions
within their corresponding local Hilbert space dimension. While it may be
theoretically straightforward to formulate and solve models with any type of
spin interactions, it is crucial to ensure that the local Hilbert space
respects the underlying Lie algebra, which physically allows for such
interactions [49, 104]. For instance, biquadratic interactions, $J_{2}$, are
absent for spin states that respect the $\mathfrak{su}(2)$ Lie algebra which
is described by Pauli matrices. To address such interactions in a physically
realistic setting, it is necessary to formulate the problem while respecting
the $\mathfrak{su}(3)$ Lie algebra, with generators in form of Gell-Mann
matrices. To illustrate the profound impact of changing the local Hilbert
space dimension and its local spin degree of freedom, our results in Fig. 1
provide a compelling demonstration using a specific example. Also other
examples, such as the BBQ model on the triangular lattice, reveals very
different properties for spins living in $\mathbb{CP}^{1}$ (classical
Heisenberg spins) [105, 106], compared to spins in $\mathbb{CP}^{2}$ [24, 39,
47]. In particular, realistic Mott insulators, in their effective low-energy
description, can accommodate hybridized multipole moments, which can be
efficiently studied within their corresponding SU($N$) representation [50,
107].

The U(3) formalism [47] has been used in the present paper as a semiclassical
approximation that neglects quantum entanglement between different sites. A
very promising and important direction for future studies would involve the
treatment of the same model with methods which can account for long-range
entanglement. In such a case, we expect that quantum fluctuations will broaden
the Kitaev SL phases, which, in our semiclassical simulations, are confined to
the singular Kitaev points located at the north and south poles of the phase
diagram in Fig. 1. In fact, this expectation has already been realized by
recent iPEPS studies in Ref. [108], which examined the unfrustrated regime of
the $S=1$ BBQ-Kitaev model for $1<\phi/\pi<2$ [cf. Fig. 1(a)]. Moreover,
recent work utilizing DMRG techniques on the $S=1$ Kitaev model has identified
the emergence of chiral ordered phases in presence of symmetric off-diagonal
interactions [109], and vortex states in the presence of single-ion anisotropy
[110]. These findings provide intriguing results that underscore the richness
of the quantum problem, offering promising avenues for further exploration.
Furthermore, given that the CSL phases in the $\mathbb{CP}^{2}$ model persist
at finite temperatures in close proximity to the Kitaev points, it might be
possible that quantum fluctuations will stabilize the CSL phases even in the
ground state. Additionally, it is known that the SO phase in the BBQ model
will transform into a plaquette valence bond solid (pVBS), as observed from
iPEPS simulations [87]. It would be interesting to investiagte how far the
pVBS extends in the presence of Kitaev interactions and whether the one-
dimensional phase (1D/q1D) survives in this scenario or not. If such a phase
persists in a quantum simulation, it could potentially stabilize exotic phases
like a Haldane phase [111], quantum loop states [112], or other phases with
dimensionally reduced characteristics, as observed in anisotropic $S=1/2$
Kitaev models [113, 114].

The $S=1$ Kitaev model with BBQ interactions offers a diverse range of
interesting and exotic phases. In the Kitaev limit we recovered the
semiclassical analog of the $S=1$ Kitaev SL, as suggested for honeycomb
materials with Ni2+ ions [71]. However, when considering the experimental
realization of higher-$S$ Kitaev materials, it becomes apparent that achieving
ideal interactions poses significant challenges. $S=1$ magnets naturally
introduce higher-order biquadratic interactions, placing potential candidate
materials likely in the vicinity of the Kitaev limit, but with finite $J_{1}$
and $J_{2}$ interactions. In such cases, the surrounding phases become
relevant for characterizing and classifying synthesized materials. A
ferroquadrupolar spin nematic state is expected for negative biquadratic
interactions, $J_{2}$, which can be expected in materials with spin-phonon
coupling [115]. States, such as triple-$q$, q1D, and CSL states appear for
positive biquadratic interactions $J_{2}$, which are expected in materials
with orbital degeneracy [116, 117, 118]. Phases with scalar spin chirality,
such as the FM/AFM CSL or triple-$q$ ordered states found in this study, are
of particular interest, as they may lead to unconventional phenomena, such as
the magnon Hall effect [119, 120, 121]. In fact, triple-$q$ states were also
found in the Kitaev-Heisenberg model under magnetic field [122] and play a
crucial role in understanding the physics of materials such as the honeycomb
cobaltites Na2Co2TeO6 [123, 124], and Na3Co2SbO6 [125]. Despite the presence
of non-Kitaev interactions, the CSL presented in our study is stable at finite
temperatures in the vicinity of the Kitaev limit. This observation suggests
the intriguing possibility of stabilizing a similar state in real materials,
potentially leading to the emergence of a novel Kitaev-like SL that breaks
time-reversal symmetry.

###### Acknowledgements.

The authors are pleased to acknowledge helpful conversations with Matthias
Gohlke, Yasuyuki Kato, Geet Rakala, and Kotaro Shimizu. This work was
supported by the Quantum Liquid Crystals JSPS KAKENHI Grants No. JP19H05822
and No. JP19H05825, MEXT as "Program for Promoting Research on the
Supercomputer Fugaku", Grant No. JPMXP1020230411, and the Theory of Quantum
Matter Unit, OIST. Numerical calculations were carried out using HPC
facilities provided by the Supercomputer Center of the Institute for Solid
State Physics, the University of Tokyo.

##  Appendix A Details of methods

In this Appendix, we describe the numerical methods used in this study.
Simulations are performed for different local Hilbert space dimensions,
namely, SU(3) spin-coherent states on the complex-projective plane
$\mathbb{CP}^{2}$, SU(2) spin-coherent states on $\mathbb{CP}^{1}$, and
discretized spins in the eight-color model. Energies in each case are
evaluated at the level of classical approximations on local bonds, and hence
do not allow us to treat long-range quantum entanglement beyond a single site.

###  A.1 Variational energy minimization

To determine the ground-state phase diagrams, as shown in Fig. 1, we perform
large-scale variational energy minimization using the machine learning library
JAX [79]. In Fig. 1(a) we minimize $\mathcal{H}^{\mathcal{A}}_{{\sf BBQ-K}}$
in Eq. (7) by optimizing $\theta_{1}$, $\theta_{2}$, $\phi_{1}$ and $\phi_{2}$
in Eq. (6) as independent variational parameters at each site on the lattice.
Similarly, in Fig. 1(b) we minimize $\mathcal{H}^{\mathcal{S}}_{{\sf BBQ-K}}$
in Eq. (1) using local variational parameters $\theta_{1}$ and $\phi_{1}$ in
Eq. (14). We employ the gradient processing and optimization library “Optax”
with the optimizer “Adam” [80]. The initial conditions considered include
relevant ordered and disordered states, as depicted in Fig. 1(d), as well as
random initial states. The optimization process involves $5000$-$10000$ steps
per model parameter.

###  A.2 Monte Carlo simulation

To obtain thermodynamic properties we perform classical MC simulations where
spins are locally updated at every site on the physical lattice of size
$N_{\rm S}$. In the $\mathbb{CP}^{2}$ spin space, we sample $\mathcal{A}$
matrices by choosing parameters $0\leq\theta_{1},\theta_{2}\leq 1$ and
$0\leq\phi_{1},\phi_{2}<2\pi$ in Eq. (6) at random, while in the
$\mathbb{CP}^{1}$ spin space we sample classical vector spins with
$-1\leq\cos{\theta_{1}}\leq 1$ and $0\leq\phi_{1}<2\pi$ in Eq. (14). For a
spin flip in the eight-color model we randomly select, with equal weight, one
out of all eight possible spin states, as defined in Table 1, Following the
standard single-spin flip Metropolis algorithm [126] we accept or reject a new
spin state at site $i$, after evaluating the local bond energies. In the
actual calculations, a single MC step consists of $N_{\rm S}$ local spin-flip
attempts on randomly chosen sites. Simulations are carried out in parallel for
replicas at different temperatures, using the replica-exchange method,
initiated by the parallel tempering algorithm [127, 128, 129] every 100 MC
steps. Results for thermodynamic quantities are averaged over $5\times 10^{5}$
statistically independent samples, after initial $5\times 10^{5}$ MC steps for
slowly heating from the ground state, as obtained by the variational energy
minimization and and further $5\times 10^{5}$ MC steps for thermalization.
Error bars are estimated by averaging over 10 independent simulation runs.

MC simulations in the frustrated parameter region, where the eight-color CSL
appears, suffer from severe slowing down, due to the failure of single-spin
flip updates. As explained in Sec III.2.1, the CSL is expected to have an
extensive number of degenerate states, which are connected via simultaneous
flips of a group of spins covering at least one hexagon (see Fig. 3). To
efficiently sample over the full manifold of states, we combine single-spin
flip updates with an additional cluster update of six spins. For the eight-
color model and the $\mathbb{CP}^{1}$ model we adopt Eq. (28), while for the
$\mathbb{CP}^{2}$ model, where spins are represented by $3\times 3$ Hermitian
$\mathcal{A}$-matrices [see Eqs. (3)–(6)], we transform Eq. (28) into the
matrix form

$\displaystyle\mathcal{A}_{p,j}^{\rm
new}=(R^{\alpha}_{j})^{-1}\mathcal{A}_{p,j}\ R^{\alpha}_{j}\,.$ (40)

By comparing energies before and after the cluster update, we accept or reject
the new state following the Metropolis argument [126]. Within one MC step we
implement $N_{\rm h}=N_{\rm S}/2$ cluster-update attempts for randomly chosen
hexagons on the whole lattice. The efficiency of this update within the CSL
phase strongly depends on the degrees of freedom of spins. As discussed for
the acceptance ratio in Appendix D, this cluster update becomes rejection free
in the eight-color model, while showing a reduced efficiency for
$\mathbb{CP}^{1}$ and $\mathbb{CP}^{2}$ model calculations.

###  A.3 Molecular dynamics simulation

Dynamics are present in models with continuous spin degree of freedom, namely
for the $\mathbb{CP}^{2}$ and $\mathbb{CP}^{1}$ models. The equations of
motion for spins respecting the spin space $\mathbb{CP}^{2}$ [spin coherent
states of SU(3)] are expressed within our formalism in terms of $\mathcal{A}$
matrices (see Sec. II). This allows us to rewrite the original BBQ-Kitaev
Hamiltonian of Eq. (1) into the bilinear form of Eq. (7),

$\mathcal{H}^{\mathcal{A}}_{{\sf BBQ-K}}=\mathcal{H}^{\mathcal{A}}_{\sf
BBQ}+\mathcal{H}^{\mathcal{A}}_{\sf K[{\it x}]}+\mathcal{H}^{\mathcal{A}}_{\sf
K[{\it y}]}+\mathcal{H}^{\mathcal{A}}_{\sf K[{\it z}]}$ (41)

with

$\displaystyle\mathcal{H}^{\mathcal{A}}_{\sf BBQ}$
$\displaystyle=\sum_{\langle
ij\rangle}\Big{[}J_{1}\mathcal{A}^{\alpha}_{i\beta}\mathcal{A}^{\beta}_{j\alpha}+(J_{2}-J_{1})\mathcal{A}^{\alpha}_{i\beta}\mathcal{A}^{\alpha}_{j\beta}+J_{2}\mathcal{A}^{\alpha}_{i\alpha}\mathcal{A}^{\beta}_{j\beta}\Big{]}\,,$
(42) $\displaystyle\mathcal{H}^{\mathcal{A}}_{\sf K[{\it x}]}$
$\displaystyle=-K_{x}\sum_{\langle
ij\rangle_{x}}\left(\mathcal{A}^{y}_{iz}-\mathcal{A}^{z}_{iy}\right)\left(\mathcal{A}^{y}_{jz}-\mathcal{A}^{z}_{jy}\right)\,,$
(43) $\displaystyle\mathcal{H}^{\mathcal{A}}_{\sf K[{\it y}]}$
$\displaystyle=-K_{y}\sum_{\langle
ij\rangle_{y}}\left(\mathcal{A}^{x}_{iz}-\mathcal{A}^{z}_{ix}\right)\left(\mathcal{A}^{x}_{jz}-\mathcal{A}^{z}_{jx}\right)\,,$
(44) $\displaystyle\mathcal{H}^{\mathcal{A}}_{\sf K[{\it z}]}$
$\displaystyle=-K_{z}\sum_{\langle
ij\rangle_{z}}\left(\mathcal{A}^{x}_{iy}-\mathcal{A}^{y}_{ix}\right)\left(\mathcal{A}^{x}_{jy}-\mathcal{A}^{y}_{jx}\right)\,,$
(45)

where $\mathcal{A}$ matrices respect the commutation relations [47]:

$[\mathcal{A}^{\alpha}_{i\beta},\mathcal{A}^{\gamma}_{i\eta}]=\delta^{\gamma}_{\beta}\mathcal{A}^{\alpha}_{i\eta}-\delta^{\alpha}_{\eta}\mathcal{A}^{\gamma}_{i\beta}\,,\quad[\mathcal{A}^{\alpha}_{i\beta},\mathcal{A}^{\gamma}_{j\eta}]=0\,.$
(46)

By using Eq. (46), one can solve the Heisenberg equations of motion, written
in terms of $\mathcal{A}$ matrices

$\frac{d}{dt}\mathcal{A}^{\alpha}_{i\beta}=\frac{d}{dt}\mathcal{A}^{\alpha}_{i\beta}\Big{|}_{\sf
BBQ}+\frac{d}{dt}\mathcal{A}^{\alpha}_{i\beta}\Big{|}_{\sf K[{\it
x}]}+\frac{d}{dt}\mathcal{A}^{\alpha}_{i\beta}\Big{|}_{\sf K[{\it
y}]}+\frac{d}{dt}\mathcal{A}^{\alpha}_{i\beta}\Big{|}_{\sf K[{\it z}]}\,,$
(47)

which can be explicitly solved for contributions from the BBQ interactions

$\displaystyle\frac{d}{dt}\mathcal{A}^{\alpha}_{i\beta}\Big{|}_{\sf BBQ}=$
$\displaystyle-i\left[\mathcal{A}^{\alpha}_{i\beta},\mathcal{H}^{\mathcal{A}}_{{\sf
BBQ}}\right]$ (48) $\displaystyle=$
$\displaystyle-i\sum_{\delta}\Big{[}J_{1}\left(\mathcal{A}^{\alpha}_{i\gamma}\mathcal{A}^{\gamma}_{i+{\delta},\beta}-\mathcal{A}^{\gamma}_{i\beta}\mathcal{A}^{\alpha}_{i+{\delta},\gamma}\right)$
$\displaystyle+(J_{2}-J_{1})\left(\mathcal{A}^{\alpha}_{i\gamma}\mathcal{A}^{\beta}_{i+{\delta},\gamma}-\mathcal{A}^{\gamma}_{i\beta}\mathcal{A}^{\gamma}_{i+{\delta},\alpha}\right)\Big{]}\,,$

and for the bond-dependent Kitaev interactions, e.g., on the $x$ bond:

$\displaystyle\frac{d}{dt}\mathcal{A}^{\alpha}_{i\beta}\Big{|}_{\sf K[{\it
x}]}=$
$\displaystyle-i\left[\mathcal{A}^{\alpha}_{i\beta},\mathcal{H}^{\mathcal{A}}_{\sf
K[{\it x}]}\right]$ (49) $\displaystyle=$ $\displaystyle\
iK_{x}\left[\mathcal{A}^{\alpha}_{i\beta},(\mathcal{A}^{y}_{iz}-\mathcal{A}^{z}_{iy})\right]\sum_{\delta}\left(\mathcal{A}^{y}_{i+{\delta},z}-\mathcal{A}^{z}_{i+{\delta},y}\right)$
$\displaystyle=$ $\displaystyle\
iK_{x}\left(\delta^{y}_{\beta}\mathcal{A}^{\alpha}_{iz}-\delta^{\alpha}_{z}\mathcal{A}^{y}_{i\beta}-\delta^{z}_{\beta}\mathcal{A}^{\alpha}_{iy}+\delta^{\alpha}_{y}\mathcal{A}^{z}_{i\beta}\right)$
$\displaystyle\times\sum_{\delta}\left(\mathcal{A}^{y}_{i+{\delta},z}-\mathcal{A}^{z}_{i+{\delta},y}\right)\,,$

where $\sum_{\delta}$ sums over the nearest neighbors of site $i$. This
explicitly gives for the $x$ bonds

$\displaystyle\begin{rcases}\frac{d}{dt}\mathcal{A}^{\alpha}_{iy}\Big{|}_{\sf
K[{\it x}]}&=\mathcal{A}^{\alpha}_{iz}\\\
\frac{d}{dt}\mathcal{A}^{\alpha}_{iz}\Big{|}_{\sf K[{\it
x}]}&=-\mathcal{A}^{\alpha}_{iy}\\\
\frac{d}{dt}\mathcal{A}^{y}_{i\beta}\Big{|}_{\sf K[{\it
x}]}&=\mathcal{A}^{z}_{i\beta}\\\
\frac{d}{dt}\mathcal{A}^{z}_{i\beta}\Big{|}_{\sf K[{\it
x}]}&=-\mathcal{A}^{y}_{i\beta}\end{rcases}\times
iK_{x}\sum_{\delta}\left(\mathcal{A}^{y}_{i+{\delta},z}-\mathcal{A}^{z}_{i+{\delta},y}\right)\,,$
(50)

while contributions on the $y$ and $z$ bonds can be obtained by cyclic
permutations.

The equations of motion for spins respecting the spin space $\mathbb{CP}^{1}$
[spin coherent states of SU(2)] are formulated in terms of classical
Heisenberg spins, see Eq. (14). We obtain from Eq. (1) the equations of motion

$\displaystyle\frac{d}{dt}S^{\alpha}_{i}=$
$\displaystyle-i[S_{i}^{\alpha},\mathcal{H}^{\mathcal{S}}_{\sf BBQ}]$ (51)
$\displaystyle=$ $\displaystyle-iJ_{1}\sum_{\langle
ij\rangle}\left[S_{i}^{\alpha},({\bf S}_{i}\cdot{\bf S_{j}})\right]$
$\displaystyle-2iJ_{2}\sum_{\langle ij\rangle}\left[S_{i}^{\alpha},\left({\bf
S}_{i}\cdot{\bf S}_{j}\right)\right]({\bf S}_{i}\cdot{\bf S}_{j})$
$\displaystyle-i\sum_{\langle
ij\rangle_{\gamma}}K_{\gamma}\left[S_{i}^{\alpha},S^{\gamma}_{i}S^{\gamma}_{j}\right]\,,$

which explicitly gives

$\frac{d}{dt}{\bf S}_{i}={\bf
S}_{i}\times\sum_{\delta}\left\\{\left[J_{1}+2J_{2}({\bf S}_{i}\cdot{\bf
S}_{i+\delta})\right]{\bf
S}_{i+\delta}+K_{\delta}S^{\delta}_{i+\delta}\right\\}\,,$ (52)

where “$\times$” corresponds to the cross product between vector spins on the
left and the local exchange field on the right. Throughout this article we
choose isotropic Kitaev interactions with $K_{x}=K_{y}=K_{z}=K$.

After sufficient thermalization at a specific temperature with MC simulations
we evolve $1000$ statistically independent $\mathcal{A}$-matrix (${\bf}S$
vector) configurations in time by numerically integrating the equations of
motion using a fourth-order Runge-Kutta (RK-4) algorithm [130, 131]. By using
a time interval of $\delta t\approx 0.03$, we obtain a time series of
$\mathcal{A}$ matrices, $\\{\mathcal{A}^{\alpha}_{i\beta}(t)\\}$, and
classical vectors, $\\{{\bf S}_{i}(t)\\}$, which conserve the total energy of
the system with a controlled error of $\mathcal{O}(\delta t^{5})$, as expected
from the RK-4 method.

##  Appendix B Measurement of observables

###  B.1 Specific heat and entropy

During the MC measurement, we compute the specific heat

$C=\frac{1}{k_{B}T^{2}}\left[\langle\mathcal{H}^{2}\rangle-\langle\mathcal{H}\rangle^{2}\right]\,,$
(53)

where $k_{B}=1$ throughout this paper, and $\langle\ldots\rangle$ represents
an average over measurements of statistically-independent replica. In the
discretized eight-color model we obtain the thermodynamic entropy, by
integrating the specific heat over the temperature

$S(T)=S(T_{\infty})-\int_{T}^{T_{\infty}}\frac{C(T^{\prime})}{T^{\prime}}dT^{\prime}\,,$
(54)

with $T_{\infty}$ the highest temperature available.

###  B.2 Correlations in momentum space

Correlations between magnetic moments in momentum space are accessible from
the equal-time structure factor

$S_{\rm{\lambda}}({\bf
q})=\left\langle\sum_{\alpha\beta}|{m_{\lambda}}^{\alpha}_{\leavevmode\nobreak\
\beta}({\bf q})|^{2}\ \right\rangle\,,$ (55)

where the index $\lambda$ denotes the channel for $\mathcal{A}$ matrices,
dipoles, or quadrupoles as $\lambda=\mathcal{A}$, ${\rm S}$, and ${\rm Q}$,
respectively. Here, ${m_{\lambda}}^{\alpha}_{\leavevmode\nobreak\ \beta}({\bf
q})$ is defined by

$\displaystyle{m_{\rm S}}^{\alpha}_{\leavevmode\nobreak\ \alpha}({\bf q})$
$\displaystyle=-i\sum_{\beta,\gamma}\epsilon^{\alpha\ \gamma}_{\
\beta}{m_{\mathcal{A}}}^{\beta}_{\leavevmode\nobreak\ \gamma}({\bf q})\,,$
(56) $\displaystyle{m_{\rm Q}}^{\alpha}_{\leavevmode\nobreak\ \beta}({\bf q})$
$\displaystyle=-{m_{\mathcal{A}}}^{\alpha}_{\leavevmode\nobreak\ \beta}({\bf
q})-{m_{\mathcal{A}}}^{\beta}_{\leavevmode\nobreak\ \alpha}({\bf
q})+\frac{2}{3}\delta^{\alpha}_{\beta}\sum_{\gamma}{m_{\mathcal{A}}}^{\gamma}_{\leavevmode\nobreak\
\gamma}({\bf q})\;,$ (57)

where the Fourier transform of the $\mathcal{A}$ matrices is given as

${m_{\mathcal{A}}}^{\alpha}_{\leavevmode\nobreak\ \beta}({\bf
q})=\frac{1}{\sqrt{N_{\rm S}}}\sum_{i}^{N_{\rm S}}e^{i{\bf q}\cdot{\bf
r}_{i}}\mathcal{\mathcal{A}}^{\alpha}_{i\beta}\,.$ (58)

For simulations in the order parameter space $\mathbb{CP}^{1}$ and the eight-
color model, we evaluate dipole correlations directly

${m_{\rm S}}^{\alpha}_{\leavevmode\nobreak\ \alpha}({\bf
q})=\frac{1}{\sqrt{N_{\rm S}}}\sum_{i}^{N_{\rm S}}e^{i{\bf q}\cdot{\bf
r}_{i}}S_{i}^{\alpha}\,,$ (59)

from spin components sampled via Eq. (14).

###  B.3 Scalar spin chirality and semiclassical $Z_{2}$ flux operator

Figure 12:  (a) Definition of the scalar spin chirality on an individual
hexagon, see Eqs. (61)–(65). (b) The semiclassical analog of the
$\mathbb{Z}_{2}$-flux operator, as defined in Eq. (66).

We define the scalar spin chirality on the honeycomb lattice with

$\kappa=\frac{2}{N_{\rm S}}\sum_{p}\kappa_{p}\,,$ (60)

where the sum is taken over all hexagons $p$. As shown in Fig. 12(a) we divide
a single hexagon into four triangles which contribute to the total scalar
chirality of the full hexagon with

$\displaystyle\kappa_{p}=\kappa_{p}^{A}+\kappa_{p}^{B}+\kappa_{p}^{C}+\kappa_{p}^{D}\,,$
(61)

where

$\displaystyle\kappa_{p}^{A}$ $\displaystyle={\bf S}_{p0}\cdot\left({\bf
S}_{p1}\times{\bf S}_{p5}\right),$ (62) $\displaystyle\kappa_{p}^{B}$
$\displaystyle={\bf S}_{p1}\cdot\left({\bf S}_{p2}\times{\bf S}_{p3}\right),$
(63) $\displaystyle\kappa_{p}^{C}$ $\displaystyle={\bf S}_{p3}\cdot\left({\bf
S}_{p4}\times{\bf S}_{p5}\right),$ (64) $\displaystyle\kappa_{p}^{D}$
$\displaystyle={\bf S}_{p1}\cdot\left({\bf S}_{p3}\times{\bf S}_{p5}\right).$
(65)

We further define the $S=1$ semiclassical analog of the
“$\mathbb{Z}_{2}$-flux” operator [75, 57]:

$W_{p}=\prod_{j\in p}\sqrt{3}S_{j}^{\alpha}$ (66)

with $\alpha$ referring to the site-dependent label for the Kitaev bond
pointing outwards the hexagon $p$ at site $j=0,\cdots,5$. A visual
representation of Eq. (66) is given in Fig. 12(b). We emphasize that $W_{p}$
differs from the conserved flux operator in the quantum $S=1$ Kitaev model
[57].

###  B.4 Dynamical structure factor

To obtain the dynamical structure factors for dipole and quadrupole channels,
we take the Fourier transform of the time-series of $\mathcal{A}$ matrices
$\\{\mathcal{A}^{\alpha}_{i\beta}(t)\\}$ (classical vectors $\\{{\bf
S}_{i}(t)\\})$, as obtained from from molecular dynamics simulations (see
Appendix A.3), from real-space into the momentum-space and frequency domain

$S_{\rm{\lambda}}({\bf
q},\omega)=\left\langle\sum_{\alpha\beta}|{m_{\lambda}}^{\alpha}_{\leavevmode\nobreak\
\beta}({\bf q},\omega)|^{2}\ \right\rangle\,,$ (67)

with $\lambda=\mathcal{A}$, ${\rm S}$, and ${\rm Q}$. The Fourier transform

${m_{\mathcal{A}}}^{\alpha}_{\leavevmode\nobreak\ \beta}({\bf
q},\omega)=\frac{1}{\sqrt{N_{t}}}\sum_{n}^{N_{t}}e^{i\omega
t_{n}}{m_{\mathcal{A}}}^{\alpha}_{\leavevmode\nobreak\ \beta}({\bf
q},t_{n})\,,$ (68)

for $\mathcal{A}$ matrices is obtained by the fast Fourier transform (FFT)
[132], after convoluting with a Gaussian envelope to avoid numerical artifacts
like the Gibbs phenomenon [133]. Equivalent equations for dipole moments,
${m_{\rm S}}^{\alpha}_{\leavevmode\nobreak\ \alpha}({\bf q},\omega)$, and
quadrupole moments, ${m_{\rm Q}}^{\alpha}_{\leavevmode\nobreak\ \beta}({\bf
q},\omega)$, are defined in analogy to Eqs. (56) and (57), respectively.

In Figs. 10, 11, and Appendix E.2, we set the total number of time steps
$N_{t}=500$ for a maximal frequency of $\omega_{\rm max}=10J$, while in
Appendix E.1, we set $\omega_{\rm max}=6J$. Data are sampled over $1000$
statistically independent replica.

The Heisenberg equations of motion [see Eqs. (47)–(LABEL:eq:EOM3) and Eq.
(52)] are derived using bosonic commutation relations [c.f. Eq. (46)],
resulting in a dispersion relation analogous to those found in flavor-wave
theory [47]. However, the molecular dynamics simulations described above are
performed for a thermal ensemble of classical states, evolving in time while
adhering to the bosonic commutation relations. This effectively leads to
dynamics of a “mixed ensemble”, encompassing classical and bosonic statistics.

To obtain a dynamical structure factor that adheres strictly to quantum
statistics, it is necessary to introduce a prefactor $\omega/T$, which allows
to “correct” the intensity distribution via

${\tilde{S}}_{\rm{\lambda}}({\bf
q},\omega)=\frac{\omega}{T}S_{\rm{\lambda}}({\bf q},\omega)\,,$ (69)

as discussed in great detail in Sec. VII of Ref. [47], and the Supplemental
Material of Ref. [134]. Throughout the present paper, all dynamical structure
factors are presented in their “corrected” form,
${\tilde{S}}_{\rm{\lambda}}({\bf q},\omega)$.

##  Appendix C Analytic solutions of phase boundaries for CSL

In this Appendix we provide an analytical derivation for the phase boundaries
of the FM CSL in the ground state of the eight-color model, and in the
$\mathbb{CP}^{1}$ model, as discussed in Sec. III.1.

###  C.1 CSL phase at singular line in the $\mathbb{CP}^{1}$ model

In the following we derive the singular line in the phase diagram where the FM
CSL becomes the ground state in the $\mathbb{CP}^{1}$ model,
$\mathcal{H}^{\mathcal{S}}_{{\sf BBQ-K}}$ in Eq. (1), see Fig. 1(b). Without
loss of generality, we consider one possible arrangement of spins in the FM
CSL ground state with nearest-neighbors (${\bf S}_{\rm g}$, ${\bf S}_{\rm r}$)
on the $x$ bond, (${\bf S}_{\rm g}$, ${\bf S}_{\rm b}$) on the $y$ bond, and
(${\bf S}_{\rm g}$, ${\bf S}_{\rm r}$) on the $z$ bond, which all respect the
bond-dependent spin constraints of Fig. 2. We hereby follow the convention of
spin states as defined in Table 1 for the eight-color states.

We first consider the energy contribution from bilinear interactions, $J_{1}$,
and biquadratic interactions, $J_{2}$. These are bond independent, which is
why we consider energies only on the $x$ bond for simplicity. By parametrizing
spins on the Bloch sphere [see Eq. (14)], with, e.g., $\theta_{\rm g}$ and
$\phi_{\rm g}$, the polar and azimuthal angles of a “green” spin, we write:

$\displaystyle{\bf S}_{\rm g}\cdot{\bf S}_{\rm r}$
$\displaystyle=\sin{\theta_{\rm g}}\sin{\theta_{\rm r}}\cos{(\phi_{\rm
g}-\phi_{\rm r})}+\cos{\theta_{\rm g}}\cos{\theta_{\rm r}}=-\frac{1}{3}\,,$
(70)

and

$\displaystyle\left({\bf S}_{\rm g}\cdot{\bf S}_{\rm r}\right)^{2}=$
$\displaystyle\ [\sin{\theta_{\rm g}}\sin{\theta_{\rm r}}\cos{(\phi_{\rm
g}-\phi_{\rm r})}]^{2}+[\cos{\theta_{\rm g}}\cos{\theta_{\rm r}}]^{2}$ (71)
$\displaystyle+\sin{\theta_{\rm g}}\sin{\theta_{\rm r}}\cos{\theta_{\rm
g}}\cos{\theta_{\rm r}}\cos{(\phi_{\rm g}-\phi_{\rm r})}$ $\displaystyle=$
$\displaystyle\ \frac{1}{9}\,.$

Considering the spin $\tilde{\bf S}_{\rm g}$, after allowing for a small
perturbation in $\theta_{\rm g}$ by $\delta\ll 1$, gives

$\displaystyle\tilde{{\bf S}}_{\rm g}\cdot{\bf S}_{\rm r}$ $\displaystyle={\bf
S}_{\rm g}\cdot{\bf S}_{\rm r}-\delta\sin{\theta_{\rm g}\cos{\theta_{\rm r}}}$
(72) $\displaystyle={\bf S}_{\rm g}\cdot{\bf S}_{\rm
r}-\frac{\sqrt{2}}{3}\delta\,.$

Equivalently we obtain:

$\displaystyle\left(\tilde{{\bf S}}_{\rm g}\cdot{\bf S}_{\rm r}\right)^{2}$
$\displaystyle=\left({\bf S}_{\rm g}\cdot{\bf S}_{\rm
r}-\delta\sin{\theta_{\rm g}}\cos{\theta_{\rm r}}\right)^{2}$ (73)
$\displaystyle=\left({\bf S}_{\rm g}\cdot{\bf S}_{\rm
r}\right)^{2}-2\left({\bf S}_{\rm g}\cdot{\bf S}_{\rm
r}\right)\delta\sin{\theta_{\rm g}}\cos{\theta_{\rm
r}}+\mathcal{O}(\delta^{2})$ $\displaystyle=\left({\bf S}_{\rm g}\cdot{\bf
S}_{\rm r}\right)^{2}+\frac{2}{3}\frac{\sqrt{2}}{3}\delta\,,$

where we neglect terms of higher order $\mathcal{O}(\delta^{2})$.

The energies from bond-dependent Kitaev interactions contribute on the $x$
bond:

$\displaystyle\tilde{{\bf S}}_{\rm g}\cdot{\bf S}_{\rm r}$ $\displaystyle={\bf
S}_{\rm g}\cdot{\bf S}_{\rm r}-\delta\cos{\theta_{\rm g}\cos{\phi_{\rm
g}}\sin{\theta_{\rm r}}\cos{\phi_{\rm r}}}$ (74) $\displaystyle={\bf S}_{\rm
g}\cdot{\bf S}_{\rm r}+\frac{\delta}{2}\frac{\sqrt{2}}{3}\,,$

on the $y$ bond:

$\displaystyle\tilde{{\bf S}}_{\rm g}\cdot{\bf S}_{\rm b}$ $\displaystyle={\bf
S}_{\rm g}\cdot{\bf S}_{\rm b}+\delta\cos{\theta_{\rm g}\sin{\phi_{\rm
g}}\sin{\theta_{\rm b}}\sin{\phi_{\rm b}}}$ (75) $\displaystyle={\bf S}_{\rm
g}\cdot{\bf S}_{\rm b}+\frac{\delta}{2}\frac{\sqrt{2}}{3}\,,$

and on the $z$ bond:

$\displaystyle\tilde{{\bf S}}_{\rm g}\cdot{\bf S}_{\rm r}$ $\displaystyle={\bf
S}_{\rm g}\cdot{\bf S}_{\rm r}-\delta\sin{\theta_{\rm g}\cos{\theta_{\rm r}}}$
(76) $\displaystyle={\bf S}_{\rm g}\cdot{\bf S}_{\rm
r}-\delta\frac{\sqrt{2}}{3}\,.$

The energy difference, $\Delta E$, between the CSL ground state and a state
with small perturbation, using Eqs. (72)–(76), is

$\displaystyle\Delta
E=\frac{\sqrt{2}}{3}\delta\left(-J_{1}+\frac{2}{3}J_{2}\right)+K\delta\left(\frac{1}{2}\frac{\sqrt{2}}{3}+\frac{1}{2}\frac{\sqrt{2}}{3}-\frac{\sqrt{2}}{3}\right)\,.$
(77)

Thus, the CSL ground state is realized when $\Delta E=0$, which occurs for

$J_{1}=\frac{2}{3}J_{2}\,.$ (78)

This condition corresponds to the model parameters, defined in Eq. (2)

$\phi/\pi=\frac{1}{\pi}\arctan{\left(\frac{3}{2}\right)}\approx 0.312833\,,$
(79)

as found from numerics in Fig. 1(b), and stated in Eq. (16). As the
contributions from Kitaev energies perfectly cancel out, the CSL phase forms a
singular straight line in the phase diagram. A similar energy comparison can
be applied to the AFM CSL phase, leading to Eq. (17).

###  C.2 Phase boundaries of CSL in the eight-color model

The phase boundaries of the FM CSL phase in the eight-color model,
$\mathcal{H}^{{\sf 8c}}_{{\sf BBQ}-K}$ in Eq. (15), can also be calculated
analytically. The FM CSL is surrounded by the AFM and AFM 3$Q$ chiral ordered
phases [see Fig. 1(c)]. Without loss of generality, we consider the local-bond
energy for only one bond in the lattice, here the $x$ bond.

In the FM CSL, by respecting the bond-dependent spin constraints in Fig. 2,
one obtains for allowed nearest-neighbors, e.g., “bright green” and “bright
red”, the local bond energy $\boldsymbol{\sigma}_{\rm
g}\cdot\boldsymbol{\sigma}_{\rm r}=-1/3$. On the other side, the AFM phase
gives for states in the ground state, e.g., “bright green” and “dark green”,
$\boldsymbol{\sigma}_{\rm g}\cdot\bar{\boldsymbol{\sigma}}_{\rm g}=-1$. To
obtain the phase boundary between AFM and FM CSL in terms of model parameters
$\phi$ and $\theta$ [see Eq. (2)], we set their energies equal and obtain the
condition

$\displaystyle-\cos{\phi}+\sin{\phi}$
$\displaystyle=-\frac{1}{3}\cos{\phi}+\frac{1}{9}\sin{\phi}\,,$ (80)

which is independent of $\theta$ and therefore results in a straight line as
seen in the phase diagram of Fig. 1(c). Solving Eq. (80) results in

$\frac{\phi}{\pi}=\frac{1}{\pi}\arctan{\left(\frac{3}{4}\right)}\leavevmode\nobreak\
\approx 0.205\,.$ (81)

In the AFM 3$Q$ chiral ordered phase, the bond energy gives, e.g., for “bright
green” and “dark blue”, $\boldsymbol{\sigma}_{\rm
g}\cdot\bar{\boldsymbol{\sigma}}_{\rm b}=1/3$. We do the same comparison for
the phase boundary between the FM CSL and the AFM 3$Q$ chiral phases and
obtain the condition

$\displaystyle\cos{\phi}\sin{\theta}$
$\displaystyle=-\cos{\phi}\sin{\theta}\,,$ (82)

which is only fulfilled for

$\cos{\phi}\sin{\theta}=0\,.$ (83)

Equation (83) is satisfied for $\theta=0$ or $\pi$, which corresponds to the
parameters the Kitaev model, and for

$\frac{\phi}{\pi}=\frac{1}{2}\,,$ (84)

where the second solution $\phi/\pi=3/2$ is physically irrelevant.

The phase boundaries described by Eqs. (81) and (84), are well visible in Fig.
1(c). The calculations for phase boundaries in the AFM CSL are
straightforward, and give the symmetric solutions from the self-duality of
$\mathcal{H}^{{\sf 8c}}_{{\sf BBQ}-K}$, as written in Eq. (19).

##  Appendix D Acceptance ratio for single-spin flip and hexagon cluster
updates

Single-spin flip MC simulations within the CSL phase suffer from severe
slowing down, making it impossible to decorrelate the MC samples using local
spin updates. To overcome this issue, we combine single-spin flip updates with
a hexagon cluster update as physically motivated from the nature of the eight-
color CSL ground state manifold (see Fig. 3). Technical details of the
implementation of the cluster update are given in Sec. III.2.1 and Appendix
A.2

Figure 13:  Temperature dependence of the acceptance ratio in the FM CSL as
obtained in (a) the eight-color model, (b) the $\mathbb{CP}^{1}$ model, and
(c) the $\mathbb{CP}^{2}$ model. Acceptance ratios for single-spin flip
updates are depicted with blue circles, while results for hexagon updates are
shown with red circles. Model parameters are the same as used in Figs. 4
(right column), 7, and 9. Simulations were performed for finite-size clusters
of $L=48$ ($N_{\rm S}=4608$).

In Figs. 13(a)–13(c) we show the temperature dependence of the acceptance
ratio for single-spin flip and hexagon cluster updates in the eight-color
model, the $\mathbb{CP}^{1}$ model, and the $\mathbb{CP}^{2}$ model,
respectively. Model parameters are chosen to stabilize the FM CSL at low
temperatures, as depicted in Figs. 4 (right column), 7, and 9. As expected
from conventional MC simulations, the efficiency of the single-spin flip
update monotonically decreases by reducing the temperature, showing a
basically vanishing acceptance ratio below the transition into the CSL. We
note that the value of $1/8$ in Fig. 13(a) comes from the fact that we accept
a new state even when we sample over the same old state.

The hexagon update follows the same trend in the high-temperature regime.
However, below the phase transition $T_{c}$ in Figs. 13(a) and 13(b) and
$T_{c1}$ in Fig. 13(c), the acceptance ratio for the hexagon update increases
and shows finite values. In fact, the hexagon update becomes rejection free in
the eight-color model, while reaching to $\approx 25\%$ in the
$\mathbb{CP}^{1}$ model. The efficiency of the cluster update in the
$\mathbb{CP}^{2}$ model is strongly suppressed ($\lesssim 1.5\%$), as shown in
the inset of Fig. 13(c), but is still efficient to decorrelate the finite-
temperature FM CSL phase.

Figure 14:  The semiclassical $S=1$ Kitaev SL in the $\mathbb{CP}^{2}$ model,
respectively at $\theta/\pi=0$ and $1$ of $\mathcal{H}^{\mathcal{A}}_{{\sf
BBQ-K}}$ in Eq. (7). The specific heat $C/N_{\rm S}$ [Eq. (53)] in (a) shows a
crossover from a high-temperature paramagnet (PM) into the semiclassical
analog of the Kitaev spin liquid (SL). The spin structure factors for dipoles,
$S_{\rm{S}}({\bf q})$, and quadrupoles, $S_{\rm{Q}}({\bf q})$, [see Eqs.
(55)–(57)] are shown respectively for the AFM Kitaev SL in (b) and (d), and
for the FM Kitaev SL in (c) and (e). Simulations were performed at $T=0.01$
[black arrow in (a)] for finite-size clusters of $L=48$ ($N=4608$).  Figure
15:  Dynamical structure factors [see Eqs. (67)–(69)] for the semiclassical
analog of the $S=1$ AFM and FM Kitaev SL in the $\mathbb{CP}^{2}$ model,
respectively at $\theta/\pi=0$ and $1$ of $\mathcal{H}^{\mathcal{A}}_{{\sf
BBQ-K}}$ in Eq. (7). (a) and (b) show dynamical structure factors of dipoles,
${\tilde{S}}_{\rm{S}}({\bf q},\omega)$, respectively in the AFM and FM Kitaev
SL states, with their energy cross sections in (d)–(g). Excitations are very
similar to known results for the $S=1/2$ Kitaev model, see Ref. [89]. (c)
shows dynamical structure factors of quadrupoles, ${\tilde{S}}_{\rm{Q}}({\bf
q},\omega)$, with energy cross sections in (h)–(i). The signals for AFM and FM
Kitaev SL are equivalent [see Fig. 14(d) and 14(e)], and therefore shown only
once. Data were obtained from molecular dynamics simulations (details in
Appendix A.3) for finite-size clusters of $L=48$ $(N=4608)$, at temperatures
$T=0.01$ [black arrow in Fig. 14(a)]. (a)–(c) are plotted along the green path
in momentum space, as indicated in the left panel of Fig. 9(d).

##  Appendix E Additional dynamics simulations

In this Appendix we provide additional results of nontrivial dynamical
structure factors for the semiclassical $S=1$ AFM and FM Kitaev SL and the FM
chiral spin liquid (CSL), as discussed in Sec. III.5,

###  E.1 Semiclassical dynamics of the $S=1$ Kitaev spin liquid

In Fig. 14, we present thermodynamic properties of the semiclassical $S=1$ AFM
and FM Kitaev models at $\theta/\pi=0$ and $\theta/\pi=1$ of
$\mathcal{H}^{\mathcal{A}}_{{\sf BBQ-K}}$ in Eq. (7). The specific heat
$C/N_{\rm S}$ in Fig. 14(a) remains identical for AFM and FM Kitaev models,
showing a crossover from a high-temperature paramagnet into a low-temperature
cooperative paramagnet. Here, and throughout the main text, we refer to this
cooperative paramagnetic phase as the semiclassical analog of the Kitaev SL.
We numerically estimate a specific heat value in the limit of $T\to 0$, as

$\frac{1}{N_{\rm S}}C(T\to 0)=1.712(2)\,,$ (85)

which is smaller than $2$, providing direct evidence of zero-energy modes in
the semiclassical ground state. The zero-energy modes are evident in the
dynamical structure factors, as illustrated in Fig. 15.

We further present the equal-time structure factors of dipoles,
$\tilde{S}_{\rm{S}}({\bf q})$, and quadrupoles, $\tilde{S}_{\rm{Q}}({\bf q})$,
respectively for the AFM Kitaev SL in Figs. 14(b) and 14(d), and for the FM
Kitaev SL in Figs. 14(c) and 14(e), at $T=0.01$. The structure factors in both
liquids are very diffuse and show in the dipole channel a strong similarity
with the $S=1/2$ Kitaev spin liquid [135, 136, 89]. Correlations for
quadrupoles show quantitatively the same structure for both the AFM and the FM
Kitaev SL states with high intensity at the Brillouin zone center and some
weak structure reaching out to the Brillouin zone corners.

In Fig. 15, we present the dynamical structure factors [see Eqs. (67)–(69)]
for the AFM and FM Kitaev SL phases at $T=0.01$, for the same model parameters
as used in Fig. 14. The structure factors for dipoles in Figs. 15(a) and 15(b)
and their energy cross sections in Figs. 15(d)–15(g) show a strong resemblance
to results for the $S=1/2$ case [89], indicating that correlations are
dominated by spin dipoles. Additionally, in Fig. 15(c), we present the
quadrupole structure factor, which is equivalent between the AFM and FM Kitaev
SL states, revealing a high-intensity zero-energy mode, consistently with the
reduced specific heat value in Eq. (85), and a flat band at $\omega/J\approx
2$. While the ground state is primarily dipolar, dipole characteristics
implicitly influence the quadrupole components and lead to dominant zero-
energy modes also in $S_{\rm{Q}}({\bf q})$. The flat band, solely associated
with quadrupole excitations, is determined by the interaction strength
$J_{2}$. Since $J_{2}=0$ in the “pure” Kitaev model, this quadrupole mode
becomes dispersionless. The energy cross section in Fig. 15(h) captures the
dominant quadrupole correlations at small $\omega$, as it reproduces the
pattern in the static structure factors of Figs. 14(d) and 14(e). A cross
section through the flat band in Fig. 15(i) reveals a perfectly constant and
featureless intensity, indicating the presence of localized excitations.

###  E.2 Dynamics of the FM CSL – comparison between $\mathbb{CP}^{2}$ and
$\mathbb{CP}^{1}$ models

(a)  FM CSL, $\mathbb{CP}^{2}$ model, $T=0.003$ (b)  FM CSL, $\mathbb{CP}^{1}$
model, $T=0.004$ (c)  FM CSL, $\mathbb{CP}^{1}$ model, $T=0.001$

(d)  (e)  (f)

Figure 16:  Dynamical structure factors for spin-dipole moments,
${\tilde{S}}_{\rm{S}}({\bf q},\omega)$ [see Eq. (69)], in the FM CSL phase.
Results in the top row are shown along the path indicated by green lines in
the left panel of Fig. 9(d), while the bottom row shows frequency-line cuts at
the $\Gamma$ and K′ momentum points. All results are presented for the same
model parameters: $\theta/\pi=0.06$ and $\phi/\pi=0.305$ [see Eq. (2)], and
are compared between $\mathbb{CP}^{2}$ model [Eqs. (7) and (47)] and
$\mathbb{CP}^{1}$ model [Eqs. (1) and (52)] calculations. (a) Results for the
$\mathbb{CP}^{2}$ model at $T=0.003$ [c.f. Fig. 10(b)]. (b) Results for the
$\mathbb{CP}^{1}$ model at a comparable temperature, $T=0.004$, exhibiting
essentially the same features as observed for the $\mathbb{CP}^{2}$ model in
(a). (c) Results in the $\mathbb{CP}^{1}$ model at very low temperatures,
$T=0.001$. Data were obtained from molecular dynamics simulations (details in
Appendix A.3) for finite-size clusters of $L=48$ $(N=4608)$.

The eight-color model describes a system where the spin degree of freedom is
discretized, and consequently lack dynamics. Fortunately, as discussed in
Secs. III.3 and III.4, the CSL is stabilized in models with continuous spin
degree of freedom, which allows us to simulate its dynamical properties (see
technical details in Appendix A.3). Here, we compare the dynamical structure
factors of the FM CSL between the $\mathbb{CP}^{2}$ model and the
$\mathbb{CP}^{1}$ model, and demonstrate that their qualitative signatures are
essentially the same.

In Fig. 16, we present the dynamical structure factors for spin-dipole
moments, ${\tilde{S}}_{\rm{S}}({\bf q},\omega)$ [see Eqs. (67) and (69)], in
the FM CSL phase for both the $\mathbb{CP}^{2}$ model and the