Multiflavor Mott insulators in quantum materials and ultracold atoms

Abstract

Mott insulators with large and active (or multiflavor) local Hilbert spaces widely occur in quantum materials and ultracold atomic systems, and are dubbed “multiflavor Mott insulators”. For these multiflavor Mott insulators, the spin-only description with the quadratic spin interactions is often insufficient to capture the major physical processes. In the situation with active orbitals, the Kugel-Khomskii superexchange model was then proposed. We briefly review this historical model and discuss the modern developments beyond the original spin-orbital context. These include and are not restricted to the 4d/5d transition metal compounds with the spin-orbit-entangled J = 3/2 quadruplets, the rare-earth magnets with two weakly-separated crystal field doublets, breathing magnets and/or the cluster and molecular magnets, et al. We explain the microscopic origin of the emergent Kugel-Khomskii physics in each realization with some emphasis on the J = 3/2 quadruplets, and refer the candidate multiflavor Mott insulators as “J = 3/2 Mott insulators”. For the ultracold atoms, we review the multiflavor Mott insulator realization with the ultracold alkaline and alkaline-earth atoms on the optical lattices. Despite a large local Hilbert space from the atomic hyperfine spin states, the system could naturally realize a large symmetry group such as the Sp(N) and SU(N) symmetries. These ultracold atomic systems lie in the large-N regime of these symmetry groups and are characterized by strong quantum fluctuations. The Kugel-Khomskii physics and the exotic quantum ground states with the “baryon-like” physics can appear in various limits. We conclude with our vision and outlook on this subject.

Introduction

There are many different ways to classify and view the physics related to correlated quantum many-body systems¹. One classification scheme may be insufficient and thus misses other complementary views, while different views could raise different types of physical questions and point to new directions of our field. Here we sketch some popular schemes and views in the field. One could classify the quantum many-body systems from the relevant and emergent phases and phase transitions. Alternatively, one could summarize various physical phenomena and classify the qualitative behaviors that may or may not be unique to particular phases, but these efforts can be quite useful experimentally and practically. One can further identify the internal structures of the underlying systems such as the intrinsic topological structures² or emergent symmetries and the related experimental signatures. One may also discuss various universal properties pertinent to certain phases or focus on the physical realization of these phases in quantum materials with the relevant physical degrees of freedom. The last view may necessarily involve a significant amount of specific physics and specific features of the degrees of freedom and the underlying quantum materials^3,4,5. The universal parts of the physics, however, are inevitably entangled with the specific physics and manifest themselves in terms of the specific degrees of freedom. The balance between universality and specifics may be strongly constrained by the specific materials instead of being determined by the more subjective purposes. Therefore, one could further classify the correlated many-body systems according to the relevant physical degrees of freedom and their interactions.

With the above thoughts, we turn to the multiflavor Mott insulators. As we have described in the abstract, these multiflavor Mott systems carry a large and active local Hilbert spaces with multiple flavors local states⁶. One traditional example of such multiflavor Mott insulating systems is the one involving both active spin and orbital degrees of freedom where the spin and orbital states play the role of flavors and the Kugel-Khomskii spin-orbital superexchange model was proposed by Kliment Kugel and Daniel Khomskii⁷. The advance that Kugel and Khomskii made beyond the Anderson’s mechanism⁸ of the superexchange spin interaction was to include the orbitals and treat them on the equal footing as the spin. Because of the complicated expression and the orbital involvement, the Kugel-Khomskii model did not receive a significant attention over the past few decades. Nevertheless, the Kugel-Khomskii physics is realistic and relevant for many physical systems^{3,4,9,10,11,12}. In the recent years, the orbital degrees of freedom and the orbital related physics are getting more attention. This is due to many factors, such as the rising of topological materials^13,14,15,16 (that often require the spin-orbit coupling to generate the topological bands), the spin-orbit-coupled correlated materials where the spin-orbit coupling is the key ingredient¹⁷, the experimental progress including the resonant X-ray scattering measurement that allows the experimental detection of the orbital structures and excitations^18,19, and so on. Therefore, it is timing to revisit the Kugel-Khomskii physics. Part of this review aims to suggest the broad applicability of the Kugel-Khomskii model in the multiflavor Mott insulating quantum materials beyond the original spin-orbital disentangled context. By “disentangled”, the spins and the orbitals are independent local variables. The new territory for the Kugel-Khomskii physics is proposed to be the multiflavor Mott insulators. This includes, but is not restricted to, the breathing materials and/or cluster and molecular magnets, the rare-earth magnets with weak crystal fields, the spin-orbital-entangled J = 3/2 local moments of transition metal compounds, etc. We will explain the realization of the multiflavor local Hilbert space in these Mott insulating systems. Beyond the quantum materials’ context, the ultracold atom systems, such as the alkaline atoms and the alkaline-earth atoms on the optical lattices, could serve as candidate systems to realize the Kugel-Khomskii physics together with their own merits with the exact high symmetry groups that are more proximate to the color-like degree of freedom in particle physics and will be discussed in the second half of the review.

The observation is that, despite the proposed multiflavor Mott insulators do not explicitly contain the orbitals as the original Kugel-Khomskii context, there exist emergent orbital-like degrees of freedom that play the role of the orbitals. For the breathing magnets and/or cluster Mott systems^{20,21,22,23,24,25,26,27,28,29,30,31}, it is the degenerate ground states of the local clusters that function as the effective orbitals. The effective Kugel-Khomskii physics appears when one considers the residual interactions between the degenerate ground states of the local clusters. These interactions lift the remaining degeneracy and create the many-body ground states. It is worthwhile to mention that, many of the correlated miore systems are actually in the cluster Mott insulating regime^32,33. For certain rare-earth magnets^{34,35,36,37,38,39}, although the orbitals are present implicitly, they are strongly entangled with the spin, and then what is meaningful is the local “J” moment. For this case, it is the degenerate energy levels that can be treated as effective orbital degrees of freedom. For the transition metal compounds with the J = 3/2 local moments that are dubbed “J = 3/2 Mott insulators”^{28,40,41,42,43,44,45}, one can group the four local states into two fictitious orbitals with one spin-1/2 moment and then naturally describe the interaction as the Kugel-Khomskii model. We carefully derive the superexchange model in terms of the effective spin and orbitals and express the model in the form of the Kugel-Khomskii interaction. The correspondence between the microscopic multipolar moments in the J = 3/2 language and the effective spin-orbital language is established. This correspondence may be useful for the mutual feedback between the understanding from different languages and views. This J = 3/2 local moments and the effective Kugel-Khomskii physics can be broadly applied to many 4d/5d Mott systems such as the Mo-base, Re-based, Os-based double perovskites⁴⁰, and can even be relevant to certain 3d transition metal compounds such as vanadates with the V⁴⁺ ions^3,10,46. We further discuss the physical consequences from these effective Kugel-Khomskii description for the multiflavor Mott insulators.

We devote the second half of the review to the ultracold atom system. The ultracold atom system has become a new frontier of condensed matter physics and provides new opportunities and platforms for exploring the emergent correlation physics. The Kugel-Khomskii physics turns out to be particularly relevant for many alkali and alkaline-earth atoms that possess large hyperfine spins and thus a large local Hilbert space^47,48,49. Like the multiflavor Mott insulating quantum materials in the previous paragraphs, the physical picture of these large hyperfine spins is very different from that of traditional magnets with large spins. Thus, we can discuss some of the common features shared by the multiflavor Mott insulating quantum materials and the ultracold alkali and alkaline-earth atomic Mott insulators. Here, the traditional magnets refer to the conventional 3d transition metal compounds without quenched orbitals³ and are used to distinguish from the multiflavor Mott insulating quantum materials. In these traditional magnets, large spins arise from the Hund’s coupling, and the spins of the localized electrons on the same lattice site are aligned to form a large spin. The leading order contribution to the couplings between different sites is the Anderson’s superexchange of a single pair of electrons, such that the fluctuation with respect of an ordered moment would be just ± 1. This is the physical origin of the 1/S-effect, in other words, as S is enlarged, the system evolves towards the classical direction. In many alkali and alkaline-earth atom systems, however, the situation is quite different. The energy scale is far below the atomic ionization energy, and exchanging a pair of fermions can completely shuffle the spin configuration among the 2S + 1 spin states. The spin states are much more delocalized in their Hilbert space. Thus, a large spin here behaves more like a large number of flavors or colors, which strongly enhances the quantum fluctuations. Since these ultracold alkali and alkaline-earth atoms naturally support high symmetries beyond the SU(2) for the traditional magnets, it is more appropriate to adopt the perspective of high symmetries (e.g., SU(N) and Sp(N) with N = 2S + 1) with the color-like degrees of freedom for these atomic states^{47,48,50,51,52,53,54}. Such high symmetries often require some fine-tuning for the multiflavor Mott insulating quantum materials except the more recent twisted multilayer graphenes and transition metal dichalcogenide heterostructures. Nevertheless, based on the spirit of Kugel-Khomskii physics, exchanging a pair of electrons for example in the “J = 3/2 Mott insulator” is also understood to be able to shuffle the effective spin and orbital states and make the system more delocalized in the Hilbert space with an enhanced quantum fluctuations⁴⁰.

This new perspective from these ultracold atom systems provides an opportunity to explore the many-body physics closely related to the Kugel-Khomskii-like physics for N = 4, and some aspects with high symmetries may even be connected to the high-energy physics. As an early progress, an exact and generic symmetry of Sp(4), or, isomorphically SO(5), was proved for the spin-3/2 alkali fermion systems^47,48. Under the fine-tuning, the Sp(4) symmetry can be augmented to SU(4). Later, the SU(N) symmetry has also been widely investigated in the alkaline-earth fermion systems due to the vanishing electron spins and the sole nuclear spin contribution^49,55,56. For instance, the alkaline-earth-like atom ¹⁷³Yb with the spin S = 5/2⁵⁴ has 6 components and thus SU(6) symmetry, and the ⁸⁷Sr atom with S = 9/2⁵³ has 10 spin components and thus SU(10) symmetry, respectively. Thus, we will review the properties of quantum magnetism with the ultracold atoms possessing the Sp(N) and SU(N) symmetries with N = 4. These are the ultracold-atom versions of the Kugel-Khomskii physics. They are characterized by various competing orders due to the strong quantum fluctuations. As an exotic example, they can exhibit the “baryon-like” physics. In an SU(4) quantum magnet, quantum spin fluctuations are dominated by the multi-site correlations, whose physics is beyond the two-site one of the SU(2) magnets as often studied in the condensed matter systems. It is exciting that in spite of the huge difference of the energy scales, the large-spin cold fermions can also exhibit similar physics to quantum chromodynamics (QCD).

The remaining parts of this review are organized as follows. We first start with a brief introduction of the traditional Kugel-Khomskii superexchange model in the multiflavored Mott insulating systems with active orbitals. We then turn the attention to explain various modern quantum materials’ realizations of the emergent Kugel-Khomskii physics and the status of the effective orbitals. After these examples, we focus the attention on the “J = 3/2 Mott insulators” and establish the emergent Kugel-Khomskii model. We then review the ultracold atoms on the optical lattices and discuss the high-symmetry models in the multiflavor Mott regimes and the emergent physics with the alkaline and alkaline earth atoms. Finally, we summarize this review.

Traditional Kugel-Khomskii exchange model

As we have remarked in the Introduction, a conventional and representative example of the multiflavor Mott insulators with a large and active local Hilbert space is the one involving active orbital degrees of freedom where the spin and orbital states play the role of the multiflavored local Hilbert space. Thus, we start the review with these spin-orbital-based Mott insulators and their superexchange interactions.

The Anderson superexchange interaction for the spin degrees of freedom in the Mott insulators is widely accepted as the major mechanism for the antiferromagnetism⁸. Anderson’s treatment was perturbative. The virtual exchange of the localized electrons from the neighboring sites through the high-order perturbation processes generates a Heisenberg interaction between the local spin moments. Anderson’s original work was based on a single-band Hubbard model where only one orbital related band is involved, and Anderson’s treatment could be well adjusted to include the orbitals from the intermediate anions. Insights from these calculations were summarized as the empirical Goodenough-Kanamori-Anderson (GKA) rules^57,58.

For the Hubbard model with multiple orbitals at the magnetic ions, the orbital necessarily becomes an active degree of freedom in addition to the spin once the system is in the Mott insulating phase with localized electrons on the lattice sites⁷. In this case, more structures are involved in the local moment formation that generates a large local Hilbert space with both spins and orbitals in the Mott regime, and the original Anderson’s spin exchange model cannot be directly applied to the multiflavor Mott insulators with the spin and orbitals here. With more active degrees of freedom in the large local physical Hilbert space, the GKA rules can no longer provide the nature and the magnitude of the exchange interactions, even the signs of the exchange interactions cannot be determined. These exchange interactions sensitively depend on the orbital configurations on each lattice site. Moreover, as the orbital is an active and dynamical degree of freedom here, the orbitals are intimately involved in the superexchange processes. Thus, in addition to the exchange of the spin quantum numbers, the perturbative superexchange processes in these multiflavor Mott insulators are able to exchange the orbital quantum numbers. In reality, both the spin and the orbital quantum numbers can be exchanged separately or simultaneously. Taking together, the full exchange Hamiltonian would involve pure spin exchange, pure orbital exchange, and the mixed spin-orbital exchange⁷. This spin-orbital exchange model is nowadays referred as “Kugel-Khomskii model”.

Closely following Kugel and Khomskii⁷, we present an illustrative derivation of the Kugel-Khomskii spin-orbital exchange model from a two-orbital Hubbard model. The Hubbard model is given^7,59,

$$\begin{array}{l}H\,=\,\mathop{\sum}\limits_{\langle ij\rangle }{t}_{ij}^{\alpha \beta }{a}_{i\alpha \sigma }^{{\dagger} }{a}_{j\beta \sigma }+\frac{U}{2}\mathop{\sum}\limits_{i}{n}_{i\alpha \sigma }{n}_{i\beta {\sigma }^{{\prime} }}(1-{\delta }_{\alpha \beta }{\delta }_{\sigma {\sigma }^{{\prime} }})\\ \qquad\,-\,\frac{1}{2}\mathop{\sum}\limits_{i,\alpha \ne \beta }{J}_{{{{\rm{H}}}}}({a}_{i\alpha \sigma }^{{\dagger} }{a}_{i\alpha {\sigma }^{{\prime} }}{a}_{i\beta {\sigma }^{{\prime} }}^{{\dagger} }{a}_{i\beta \sigma }+{a}_{i\alpha \sigma }^{{\dagger} }{a}_{i\beta \sigma }{a}_{i\alpha {\sigma }^{{\prime} }}^{{\dagger} }{a}_{i\beta {\sigma }^{{\prime} }}),\end{array}$$

(1)

where α, β = 1, 2 label the two orbitals, $\sigma ,{\sigma }^{{\prime} }$ label the spin quantum number, and the final term takes care of the inter-orbital pair hopping and the Hund’s coupling (J_H). The orbital degeneracy is assumed, and an isotropic and diagonal hopping t₁₁ = t₂₂ = t, t₁₂ = 0 is further assumed. In reality, the isotropic and diagonal hoppings are not guaranteed, and the orbital degeneracy is not quite necessary. A standard perturbation treatment yields the conventional Kugel-Khomskii model with

$${H}_{{{{\rm{KK}}}}}=\mathop{\sum}\limits_{\langle ij\rangle }{J}_{1}\,{{{{\boldsymbol{S}}}}}_{i}\cdot {{{{\boldsymbol{S}}}}}_{j}+{J}_{2}\,{{{{\boldsymbol{\tau }}}}}_{i}\cdot {{{{\boldsymbol{\tau }}}}}_{j}+4{J}_{3}\,({{{{\boldsymbol{S}}}}}_{i}\cdot {{{{\boldsymbol{S}}}}}_{j})\,({{{{\boldsymbol{\tau }}}}}_{i}\cdot {{{{\boldsymbol{\tau }}}}}_{j}),$$

(2)

where we have the relations between the electron billinear operators and the spin-orbital operators,

$$\mathop{\sum}\limits_{\sigma }{a}_{i,1,\sigma }^{{\dagger} }{a}_{i,1,\sigma }=\frac{1}{2}+{\tau }_{i}^{z},$$

(3)

$$\mathop{\sum}\limits_{\sigma }{a}_{i,2,\sigma }^{{\dagger} }{a}_{i,2,\sigma }=\frac{1}{2}-{\tau }_{i}^{z},$$

(4)

$$\mathop{\sum}\limits_{\sigma }{a}_{i,1,\sigma }^{{\dagger} }{a}_{i,2,\sigma }={\tau }_{i}^{+},$$

(5)

$$\mathop{\sum}\limits_{\sigma }{a}_{i,2,\sigma }^{{\dagger} }{a}_{i,1,\sigma }={\tau }_{i}^{-},$$

(6)

and

$$\mathop{\sum}\limits_{\alpha }{a}_{i,\alpha ,\uparrow }^{{\dagger} }{a}_{i,\alpha ,\uparrow }=\frac{1}{2}+{S}_{i}^{z},$$

(7)

$$\mathop{\sum}\limits_{\alpha }{a}_{i,\alpha ,\downarrow }^{{\dagger} }{a}_{i,\alpha ,\downarrow }=\frac{1}{2}-{S}_{i}^{z},$$

(8)

$$\mathop{\sum}\limits_{\alpha }{a}_{i,\alpha ,\uparrow }^{{\dagger} }{a}_{i,\alpha ,\downarrow }={S}_{i}^{+},$$

(9)

$$\mathop{\sum}\limits_{\alpha }{a}_{i,\alpha ,\downarrow }^{{\dagger} }{a}_{i,\alpha ,\uparrow }={S}_{i}^{-},$$

(10)

and the exchange couplings are given as

$${J}_{1}=\frac{2{t}^{2}}{U}\left(1-\frac{{J}_{\,{{\mbox{H}}}\,}^{}}{U}\right),$$

(11)

$${J}_{2}={J}_{3}=\frac{2{t}^{2}}{U}\left(1+\frac{{J}_{\,{{\mbox{H}}}\,}^{}}{U}\right).$$

(12)

Here S_i and τ_i refer to the electron spin and the orbital pseudospin, respectively. In particular, when τ^z = 1/2 (−1/2), the orbital 1 (2) is occupied. Because the spin and orbitals are disentangled, the spin sector retains the SU(2) rotational symmetry. The orbital SU(2) symmetry in Eq. (2) is accidental and is due to the special choice of the hopping parameters. Via a fune-tuning of the hoppings and the interactions, the Kugel-Khomskii model could have a larger symmetry such as SU(4)⁶⁰, and the limit with a higher symmetry may provide a new solvability of this complicated but realistic model.

In the above derivation of the Kugel-Khomskii model, the orbital degeneracy is not actually required. As long as the crystal field splitting between the orbitals is small or comparable to the exchange energy scale, one needs to seriously include the orbitals into the description of the correlation physics for the local moments in the Mott regime. This would clear up the unnecessary constraint of the Kugel-Khomskii physics for the Mott systems to have an explicit orbital degeneracy. A perfect orbital degeneracy requires a high crystal field symmetry and is not quite common. The clearance of this constraint would already expand the applicability of the Kugel-Khomskii physics in the transitional metal compounds.

From the above example, one could readily extract some of the essential properties for the Kugel-Khomskii model and the related physics for the multiflavor Mott insulators with active spins and orbitals. As the orbitals have the spatial orientations in the real space, the electron hopping between the orbitals from the neighboring sites depends strongly on the bond orientation and the orbital configurations. The resulting Kugel-Khomskii model is anisotropic in the orbital sector, and the exchange interaction depends on the bond orientations⁷. Because the spin and the orbital are disentangled in the parent Hubbard model and the spin-orbit coupling is not considered, the spin interaction in the spin sector, however, remains isotropic in the Kugel-Khomskii model and is described by the conventional Heisenberg interaction. For the relevant transition metal compounds, the common orbital degrees of freedom can be e_g and t_2g orbitals. In the case of the e_g orbitals, an effective pseudospin-1/2 operator is often used to label the orbital state⁷, and this corresponds to the situation in the above example. For the t_2g orbitals, a pseudospin-1 operator is used to operate on the three t_2g orbitals^61,62, and the corresponding Kugel-Khomskii model is much more involved than the e_g case. The local physical Hilbert space is significantly enlarged by the spin and orbital states. The Kugel-Khomskii model operates quite effectively in this enlarged local Hilbert space. As it was mentioned more generally in the Introduction, this operation is quite different from the spin-only Mott insulators with a large-S local spin moment where the orbital degree of freedom is quenched. Although the physical Hilbert space is enlarged with a large-S local spin moment, the simple Heisenberg spin model merely changes the spin quantum number by 1 with one operation. In contrast, the Kugel-Khomskii spin-orbital model is more effectively in delocalizing the spin-orbital states in the enlarged physical Hilbert space and thus enhances the quantum fluctuations. The fluctuation-driven orders and/or exotic quantum states can be stabilized in this fashion^5,12,63.

Modern quantum material realization of emergent Kugel-Khomskii physics

In the previous part, we have demonstrated that the traditional Kugel-Khomskii model can be derived from an extended Hubbard model with multiple orbitals for the multiflavor Mott insulators. The orbital degeneracy is often assumed in the literature but not necessarily required. If the orbital separation is of the similar energy scale as the superexchange energy scale, then one should carefully include these orbitals into the modeling. Depending on the electron filling, the spin moment can be spin-1/2 and spin-1, and the pseudospin for the orbital sector can be pseudospin-1/2 (for e_g degeneracy, or two-fold t_2g degeneracy) and pseudospin-1 moment (for three-fold t_2g degeneracy). The spin and orbitals are disentangled in this model.

A well-known example would be the Fe-based superconductors^64,65,66. Although the parent materials behave mostly like a bad metal and are thus modeled by an extended Hubbard model with multiple electron orbitals, many important physics such as the magnetic excitations and spectra may be better understood from the local moment picture. This was used to interpret the magnetic excitations in FeSe that is believed to be the most “Mott”-like Fe-based superconducting system^67,68. The active orbital degree of freedom in FeSe, however, has not been included into the analysis of the magnetic properties and the excitations⁶⁹. Thus, FeSe can be a good application of the original Kugel-Khomskii model for the understanding of the magnetism, the orbital physics and the nematicity^70,71,72.

While the Kugel-Khomskii model is proposed for the multiflavor Mott insulators with the disentangled spin and orbitals, we explain its broad application to other systems below, where the Kugel-Khomskii physics emerges from the multiflavor Hilbert space with effective orbitals and thus differs from the traditional Kugel-Khomskii model.

Breathing magnets and cluster magnets

Breathing magnets (and cluster magnets) represent a new family of magnetic materials whose building blocks are not the magnetic ions, and can find their applications in many organic magnets^73,74 and even inorganic compounds⁷⁵. Instead, the systems consist of the magnetic cluster units as the building blocks, and these magnetic cluster units provide the elementary and local degrees of freedom for the magnetism. To understand the physics of these systems, one ought to first understand the local physics on the cluster unit and find the relevant low-energy states. From the perspective of the multiflavor Mott insulators, the local states of the magnetic cluster unit contribute to the multiple flavors of states in the local Hilbert space. The many-body model for the system should be constructed from these relevant local low-energy states. To illustrate the point above, we notice that the early spin liquid candidate κ-(ET)₂Cu₂(CN)₃ can actually be placed into the category of the cluster magnets⁷⁶. In κ-(ET)₂Cu₂(CN)₃, each (ET)₂ molecular dimer hosts odd number of electrons. As the molecular dimers form a triangular lattice, it was proposed that this model realizes the triangular lattice Hubbard model at the half filling. The basis of the Hamiltonian is the Wannier functions associated with the antibonding states of the highest occupied molecular orbitals on each (ET)₂ dimer⁷⁷. More generally, the local energy levels of the magnetic clusters should be understood or classified from the irreducible representation of the local symmetry group, and the effective orbital degrees of the freedom on the cluster is then interpreted as the local basis of the irreducible representation^20,21,28.

To deliver the idea of the effective orbital degree of freedom, we start from the breathing kagomé lattice (see Fig. 1), and assume the simple Heisenberg interactions with alternating couplings,

$$H=J\mathop{\sum}\limits_{\langle ij\rangle \in \bigtriangleup }{{{{\boldsymbol{S}}}}}_{i}\cdot {{{{\boldsymbol{S}}}}}_{j}+{J}^{{\prime} }\mathop{\sum}\limits_{\langle ij\rangle \in \bigtriangledown }{{{{\boldsymbol{S}}}}}_{i}\cdot {{{{\boldsymbol{S}}}}}_{j},$$

(13)

where “△” ("▽”) refers to the up (down) triangles. The word “breathing” refers to the fact that the “△” triangles are of different size from the “▽” triangles and was first used in the context of the breathing pyrochlore magnets⁷⁸. In the strong breathing limit with $J\gg {J}^{{\prime} }$, one should first consider the local states on the up triangles and couple these local states together through the ${J}^{{\prime} }$-links on the down triangles. On each up triangles, there are three spin-1/2 local moments. With the antiferromagnetic interactions, the ground states have four fold degeneracies. This can be understood simply from the spin multiplication relation with

$$\frac{1}{2}\otimes \frac{1}{2}\otimes \frac{1}{2}=\frac{1}{2}\oplus \frac{1}{2}\oplus \frac{3}{2},$$

(14)

where the left hand side refers to the three spin local moments on the up triangle and the right hand side refers to the spin quantum number of the total spin of the up triangular cluster. The total spin, S_tot = 3/2, is only favored if the interaction is ferromagnetic. The antiferromagnetic spin interaction favors a total spin S_tot = 1/2, that is realized by forming a spin singlet between two spins and leaving the remaining spin as a dangling spin (see Fig. 1). This would naively lead to three singlet occupation configurations. It turns out that only two of them are linearly independent. Counting the spin-up and spin-down degeneracy of the dangling spin, there are in total four fold ground state degeneracies on each up triangular cluster. To make connection with the Kugel-Khomskii physics, one can simply regard the total spin of the up-triangle cluster as the effective spin, and regard the two-fold degeneracy of the spin singlet configuration as the effective orbital. The later has the spirit of an orbital as both are even under time reversal. The effective spin-orbital states on the up-triangle cluster are listed as follows,

$$\begin{array}{l}\left\vert {s}^{z}=\uparrow ,{\tau }^{z}=\uparrow \right\rangle =\frac{1}{\sqrt{3}}\left[\left\vert {\downarrow }_{1}{\uparrow }_{2}{\uparrow }_{3}\right\rangle +{e}^{i\frac{2\pi }{3}}\left\vert {\uparrow }_{1}{\downarrow }_{2}{\uparrow }_{3}\right\rangle +{e}^{-i\frac{2\pi }{3}}\left\vert {\uparrow }_{1}{\uparrow }_{2}{\downarrow }_{3}\right\rangle \right],\\ \left\vert {s}^{z}=\uparrow ,{\tau }^{z}=\downarrow \right\rangle =\frac{1}{\sqrt{3}}\left[\left\vert {\downarrow }_{1}{\uparrow }_{2}{\uparrow }_{3}\right\rangle +{e}^{-i\frac{2\pi }{3}}\left\vert {\uparrow }_{1}{\downarrow }_{2}{\uparrow }_{3}\right\rangle +{e}^{i\frac{2\pi }{3}}\left\vert {\uparrow }_{1}{\uparrow }_{2}{\downarrow }_{3}\right\rangle \right],\\ \left\vert {s}^{z}=\downarrow ,{\tau }^{z}=\uparrow \right\rangle =\frac{1}{\sqrt{3}}\left[\left\vert {\uparrow }_{1}{\downarrow }_{2}{\downarrow }_{3}\right\rangle +{e}^{-i\frac{2\pi }{3}}\left\vert {\downarrow }_{1}{\uparrow }_{2}{\downarrow }_{3}\right\rangle +{e}^{i\frac{2\pi }{3}}\left\vert {\downarrow }_{1}{\downarrow }_{2}{\uparrow }_{3}\right\rangle \right],\\ \left\vert {s}^{z}=\downarrow ,{\tau }^{z}=\downarrow \right\rangle =\frac{1}{\sqrt{3}}\left[\left\vert {\downarrow }_{1}{\uparrow }_{2}{\uparrow }_{3}\right\rangle +{e}^{-i\frac{2\pi }{3}}\left\vert {\uparrow }_{1}{\downarrow }_{2}{\uparrow }_{3}\right\rangle +{e}^{i\frac{2\pi }{3}}\left\vert {\uparrow }_{1}{\uparrow }_{2}{\downarrow }_{3}\right\rangle \right].\end{array}$$

The above four states define the four flavors of the local states for the local triangular unit if one views the system as the multiflavor Mott insulator with the up-triangle cluster as the effective lattice site. One then includes the ${J}^{{\prime} }$ interaction between the up-triangular clusters, and the resulting model is a Kugel-Khomskii model that is of the following form⁷⁹,

$$\begin{array}{l}{l}_{{{{\rm{KK}}}}}\,=\;\frac{4{J}^{{\prime} }}{9}\mathop{\sum}\limits_{\langle {{{\boldsymbol{r}}}}{{{{\boldsymbol{r}}}}}^{{\prime} }\rangle }({{{{\boldsymbol{s}}}}}_{{{{\boldsymbol{r}}}}}\cdot {{{{\boldsymbol{s}}}}}_{{{{{\boldsymbol{r}}}}}^{{\prime} }})\left[\frac{1}{2}-({\alpha }_{{{{\boldsymbol{r}}}}{{{{\boldsymbol{r}}}}}^{{\prime} }}{\tau }_{{{{\boldsymbol{r}}}}}^{-}+{\alpha }_{{{{\boldsymbol{r}}}}{{{{\boldsymbol{r}}}}}^{{\prime} }}^{* }{\tau }_{{{{\boldsymbol{r}}}}}^{+})\right]\\ \quad\qquad\times \left[\frac{1}{2}-({\beta }_{{{{\boldsymbol{r}}}}{{{{\boldsymbol{r}}}}}^{{\prime} }}^{}{\tau }_{{{{{\boldsymbol{r}}}}}^{{\prime} }}^{-}+{\beta }_{{{{\boldsymbol{r}}}}{{{{\boldsymbol{r}}}}}^{{\prime} }}^{* }{\tau }_{{{{{\boldsymbol{r}}}}}^{{\prime} }}^{+})\right],\end{array}$$

(15)

where r refers to the center of the up-triangular cluster, and s_r defines the total spin on the up-triangular cluster at r. The parameters ${\alpha }_{{{{\boldsymbol{r}}}}{{{{\boldsymbol{r}}}}}^{{\prime} }}$ and ${\beta }_{{{{\boldsymbol{r}}}}{{{{\boldsymbol{r}}}}}^{{\prime} }}$ are the bond-dependent phase factors that are consistent with the orbital-like nature of the pseudospin τ. It is found that⁷⁹, the factor ${\alpha }_{{{{\boldsymbol{r}}}}{{{{\boldsymbol{r}}}}}^{{\prime} }}$ equals to 1, e^i4π/3, or e^i2π/3 when the ${J}^{{\prime} }$-coupled bond connects two neighboring up-triangles at r and ${{{{\boldsymbol{r}}}}}^{{\prime} }$ from the spin 1, 2, 3 on the r triangle, respectively. Similarly, the factor ${\beta }_{{{{\boldsymbol{r}}}}{{{{\boldsymbol{r}}}}}^{{\prime} }}$ equals to 1, e^i4π/3, or e^i2π/3 when the ${J}^{{\prime} }$-coupled bond connects two neighboring up-triangles at r and ${{{{\boldsymbol{r}}}}}^{{\prime} }$ to the spin 1, 2, 3 on the ${{{{\boldsymbol{r}}}}}^{{\prime} }$ triangle, respectively. Since the centers of the up-triangular clusters form a triangular lattice, the emergent Kugel-Khomskii model is then defined on the triangular lattice.

Fig. 1: The breathing kagomé lattice and the spin singlet configurations.

a The breathing kagomé lattice structure with an alternating exchange coupling on the triangular cluster. b The three configurations of the spin singlet occupation on the triangular cluster. The (blue) dimer refers to the spin singlet of two spins, and the uncovered site is the dangling spin moment. Two of the three configurations are linearly independent. There are totally four ground states of the triangular cluster after including the two fold degeneracy of the dangling spin.

Apart from the breathing kagomé magnet here, the breathing pyrochlore magnet can also be understood in a similar fashion. A recent interest of the breathing pyrochlore magnet is Ba₃Yb₂Zn₅O₁₁ that is in the strong breathing limit. Due to the strong spin-orbit coupling of the Yb 4f electrons, the local Hamiltonian on the smaller tetrahedron is not a simple Heisenberg model. Nevertheless, the understanding from the irreducible representation of the tetrahedral group should still be applicable and has already been applied to the experiments³¹. The similar breathing structure occurs in the lacunar spinels^80,81 where many distinct and interesting properties emerge from the clusterization of electrons and the effective J-moments on the tetrahedra.

Rare-earth magnets with weak crystal field

It is a bit hard to imagine that the rare-earth magnets are described by the Kugel-Khomskii model. Usually, the rare-earth local moments are described by some effective spin-1/2 degrees of freedom, and this two-fold degeneracy is protected by time reversal symmetry and Kramers theorem for the Kramers doublet, and by the point group symmetry for the non-Kramers doublet. For the rare-earth local moments, the orbital degrees of freedom have already been considered from the atomic spin-orbit coupling that entangles the atomic spin with the orbitals and leads to total moment “J”. The effective spin-1/2 doublet arises from the crystal field ground state levels of the total moment J. The low-energy magnetic physics is often understood from the interaction between the effective spin-1/2 local moments. This paradigm works rather well for the pyrochlore rare-earth magnets and the triangular lattice rare-earth magnets^{82,83,84,85,86}. The reason for the success of the paradigm is due to the large crystal field gap between the ground state doublet and the excited ones in the relevant materials. If this precondition breaks down, then we need to think about other resolution.

For the Tb³⁺ ion in Tb₂Ti₂O₇ and Tb₂Sn₂O₇^{34,35,36,37,38,87,88,89}, it is known that the crystal field energy gap between the ground state doublet and the first excited doublet is not very large compared to the Curie-Weiss temperatures in these systems (see Fig. 2). Thus the ground state doublet description is insufficient to capture the low temperature magnetic properties. This regime is quoted as “weak crystal field magnetism” in ref. ³⁹. Both the ground state doublet and the first excited doublet in the left energy diagram of Fig. 2 should be included into the microscopic model, and these low-energy doublets contribute to the multiple flavors of local Hilbert space in this rare-earth version of the multiflavor Mott insulator. To think along the line of the Kugel-Khomskii physics, we assign the energy levels with the effective spin and the effective orbital configurations according to Fig. 2. Here the two effective orbitals are separated by a crystal field energy gap. The exchange interaction between the local moments would be of the Kugel-Khomskii-like. We expect other rare-earth magnets beyond Tb₂Ti₂O₇ and Tb₂Sn₂O₇ could share a similar physics from the perspective of Kugel-Khomskii and the multiflavor Mott insulators.

Fig. 2: The assignment of the effective spin and orbital to the local crystal field states of rare-earth moments.

This applies to the Tb³⁺ ion in Tb₂Ti₂O₇, Tb₂Sn₂O₇ and others. The left depicts the crystal field energy levels of the Tb³⁺ ion, and the right assigns the effective spin and orbital states.

J = 3/2 Mott insulator

What is “J = 3/2 Mott insulator”? To present this notion, we begin with the notion of “J = 1/2 Mott insulator” that seems to be quite popular in recent years^62,90,91,92. The J = 1/2 Mott insulator was proposed to be relevant to various iridates, α-RuCl₃⁹³, and even the Co-based 3d transition metal compounds^{94,95,96,97,98}. This can be understood from the Ir⁴⁺ ion under the octahedral crystal field environment^62,91. The t_2g and e_g levels for single electron states are splited by a large crystal field gap. When the spin-orbit coupling (SOC) is switched on, the t_2g orbital is entangled with the spin degree of freedom, leading to an upper J = 1/2 doublet and a lower J = 3/2 quadruplet (see Fig. 3). The Ir⁴⁺ ion has a 5d⁵ electronic configuration such that the lower quadruplet is fully filled and the upper doublet is half-filled. In the Mott insulating phase, the local moment is simply described by the spin-orbit-entangled J = 1/2 doublet. The candidate materials are often referred as “J = 1/2 Mott insulators”. As the orbital is implicitly involved into the moment, the exchange interaction inherits the orbital character and depends on the bond orientations and the components of the moments. Thus, the exchange interaction is usually not of Heisenberg like. A consequence of this anisotropic interaction is the Kitaev interaction that was popular and led to the development of the field of Kitaev materials. The major advantage of “J = 1/2 Mott insulators” on the model side is to provide extra interactions beyond the simple Heisenberg interaction, and these extra interactions could provide more opportunities to realize interesting quantum phases and orders. While it may be fashionable to name these extra interactions as Kitaev interactions and/or others, we stick to the old convention by Moriya so that the comparison can be a bit more insightful. The exchange interaction for the spin-1/2 operators is always pairwise and quadratic. According to Moriya^3,99, these interactions are classified as (symmetric) Heisenberg interaction, (antisymmetric) Dzyaloshinskii-Moriya interaction, and (symmetric) pseudo-dipole interaction.

Fig. 3: The splitting of the single electron state under the spin-orbit coupling.

The left is the crystal field scheme in the octahedral environment with the cubic symmetry. The right is the crystal field scheme under the inclusion of the spin-orbit coupling (SOC). The energy separation between the J = 1/2 doublet and the J = 3/2 quadruplet is set by the SOC.

The “J = 3/2 Mott insulator” is realized when one single electron is or three electrons are placed on the lower t_2g quadruplet and the system becomes Mott insulating^{28,40,41,42,43,44,45} as well as the lacunar spinels in the cluster localization limit⁸⁰. Despite the popularity of the “J = 1/2 Mott insulators”, the “J = 3/2 Mott insulators” did not receive much attention. We do not have any bias towards either of them, and simply address and explain what the nature could provide to us. What are the new features of the “J = 3/2 Mott insulators” from the model side? First of all, J = 3/2 local moments provide a larger local Hilbert space with more flavors of local states than J = 1/2 and allow more possibilities for interesting and unusual quantum phases and orders, thereby making the “J = 3/2 Mott insulator” an perfect example of the multiflavor Mott insulator. It was conventionally believed that large spin local moments like the spin 3/2 tend to behave more classically. This conventional belief, however, does not really apply to J = 3/2 Mott insulators. More specifically, this would come to our second point below. The conventional Heisenberg model does not generate strong quantum fluctuations as it only changes the spin quantum number by ± 1, and thus the large spin magnets with a simple Heisenberg model usually behave rather classically. For the J = 3/2 Mott insulators, more operators beyond J^x, J^y, J^z are generated in the superexchange processes and interactions due to the inclusion of the orbital degrees of freedom in the J = 3/2 local moments via the SOC. These operators are actually generators of the SU(4) group that are “4 × 4” Γ matrices. These Γ matrix operators allow the system to fluctuate more effectively between different spin states and generate stronger quantum fluctuations. This point is similar to the one that has been invoked more generally in the introduction part. Thirdly, similar to the “J = 1/2 Mott insulators”, the exchange interaction of the “J = 3/2 Mott insulators” is highly anisotropic and depends on the bond orientation.

Finally, we make a remark that the Γ matrix model for the “J = 3/2 Mott insulators” can be regarded as a parent model for various anisotropic models for the effective spin-1/2 local moments. This is understood by introducing the single-ion anisotropic term and reducing the Γ matrix model into the low-energy manifold favored by the single-ion spin term. This process may be made an analogy with the models for the free-electron band structure topology. The Luttinger model^{100,101,102,103,104,105,106,107}, that uses the Γ matrices for the k. p theory and gives rise to the Luttinger semimetal with a quadratic band touching in three dimensions, can be regarded as a parent model for generating other models for three-dimensional (3D) topological insulators and 3D Weyl semimetal upon introducing strains and magnetism.

In the next part, we will focus on the J = 3/2 Mott insulator on the triangular lattice structure and explicitly derive the Kugel-Khomskii model for this multiflavor Mott insulator.

Kugel-Khomskii model in J = 3/2 Mott insulator

Even though the J = 3/2 Mott insulator widely exists in many materials, we here consider a J = 3/2 Mott insulator in the triangular lattice for the specific demonstration purpose. It is relevant for the [111] interface of the transition metal oxide heterostructure. As we have explained in the previous part, the requirement for the magnetic ions is to have a 4d¹ or 5d¹ electron configuration in an octahedral environment where the SOC is active. To build up a physical model for this triangular lattice J = 3/2 Mott insulator, we first specify the degree of freedom. Before including the effect of the SOC, the local moment is described by a localized spin-1/2 electron on the three-fold degenerate t_2g orbitals at each magnetic ion site. The interaction between the local moments is understood as a standard Kugel-Khomskii superexchange model for the t_2g systems. Including the atomic SOC, we express the full model as

$$H={H}_{{{{\rm{KK}}}}}+{H}_{{{{\rm{SOC}}}}},$$

(16)

where the first term describes the standard Kugel-Khomskii superexchange interaction, and the second term is the on-site atomic SOC. For the z bond in Fig. 4(a), the superexchange interaction has the following form,

$${H}_{{{{\rm{KK}}}}}^{z}=\mathop{\sum}\limits_{\langle ij\rangle }J\,\left({{{{\boldsymbol{S}}}}}_{i,xy}\cdot {{{{\boldsymbol{S}}}}}_{j,xy}+\frac{1}{4}{n}_{i,xy}{n}_{j,xy}\right),$$

(17)

where S_i,xy defines the electron spin on the xy orbital with S_i,xy = S_i n_i,xy, and n_i,xy defines the electron occupation number on the xy orbital. As the xy orbital gives a dominant σ-bonding for the z bond, the superexchange interaction is primarily given by the exchange process along this σ-bonding. The electron hoppings to other orbitals or between other orbitals are expected to be quite weak and are neglected here. To keep things simple, the other Kanamori parameters such as the Hund’s coupling and the pair hopping are further neglected as they are comparably weaker than the Hubbard interaction. These approximations result in the primary contribution in the above equation for our demonstration purpose. The exchange interactions on the x and y bonds can be written down from a simple cubic permutation. Interestingly, the Kugel-Khomskii superexchange interaction for the disentangled spin-orbital states in Eq. (17) is not quite related to the emergent Kugel-Khomskii physics for the spin-orbit-entangled J = 3/2 quadruplets that is shown below.

Fig. 4: The interfacial triangular lattice of the [111] transition metal oxide heterostructure.

a The triangular lattice with “x, y, z” bond assignment for the nearest neighbors. b The [111] layer of the transition metal oxide interfaces. The coordinate system defines the coordinate system for the spin and the orbitals.

The atomic SOC has the expression,

$${H}_{{{{\rm{SOC}}}}}=\mathop{\sum}\limits_{i}-\lambda \,{{{{\boldsymbol{l}}}}}_{i}\cdot {{{{\boldsymbol{S}}}}}_{i},$$

(18)

where l_i is the effective orbital angular momentum for the t_2g orbital with l = 1, and S_i is the spin-1/2 operator for the localized electron. As it is shown in Fig. 3, the H_SOC term splits the 6-fold spin and orbital degeneracy of the single-electron states into the upper doublet with J = 1/2 and the lower quadruplets with J = 3/2. To demonstrate how the Kugel-Khomskii physics emerges in this context, we need to compare the energy scales of SO coupling λ and the superexchange J. In many 4d/5d systems, the atomic SOC is often of similar energy scales as the electron bandwidth and the electron correlation. In the Mott insulating regime, however, what is available for a meaningful comparison is the superexchange coupling. This coupling is usually much weaker than the atomic SOC for the 4d/5d materials. In the sense of perturbation theory, the atomic SOC is treated as the main Hamiltonian and the superexchange is regarded as the perturbative term. For the J = 3/2 Mott insulator with one electron per site, the atomic SOC entangles the orbital angular momentum l with the spin S, leading to a J = 3/2 quadruplet for each site. The superexchange interaction then operates on the degenerate manifold of the J = 3/2 local moments.

Following the spirit of the degenerate perturbation theory, we project the superexchange model on the degenerate manifold,

$${H}_{{{{\rm{eff}}}}}=\mathop{\prod}\limits_{i}\mathop{\sum}\limits_{{m}_{i}}\left\vert {m}_{i}\right\rangle \left\langle {m}_{i}\right\vert \cdot {H}_{{{{\rm{KK}}}}}\cdot \mathop{\prod}\limits_{j}\mathop{\sum}\limits_{{m}_{j}}\left\vert {m}_{j}\right\rangle \left\langle {m}_{j}\right\vert ,$$

(19)

where m_i (m_j) is the quantum number of the ${J}_{i}^{z}$ (${J}_{j}^{z}$) operator and takes the values of ± 1/2, ± 3/2. As we explained in the previous part, the effective Hamiltonian can be expressed into a quadratic form in terms of the Γ matrices at each site. The Γ-matrix expression, however, hides the original physical meaning. In the following, we first express the effective model in the J-basis, and then explain the emergent Kugel-Khomskii physics by expressing it into the form of Kugel-Khomskii interaction. In terms of the J operators, the effective model has the expression,

$${H}_{{{{\rm{eff}}}}}^{z}=\mathop{\sum}\limits_{\langle ij\rangle }J\,\left({\tilde{{{{\boldsymbol{S}}}}}}_{i,xy}\cdot {\tilde{{{{\boldsymbol{S}}}}}}_{j,xy}+\frac{1}{4}{\tilde{n}}_{i,xy}\cdot {\tilde{n}}_{j,xy}\right),$$

(20)

where we have

$${\tilde{S}}_{i,xy}^{x}={P}_{\frac{3}{2}}\,{S}_{i,xy}^{x}\,{P}_{\frac{3}{2}}=\frac{{J}_{i}^{x}}{4}-\frac{{J}_{i}^{z}{J}_{i}^{x}{J}_{i}^{z}}{3},$$

(21)

$${\tilde{S}}_{i,xy}^{y}={P}_{\frac{3}{2}}\,{S}_{i,xy}^{y}\,{P}_{\frac{3}{2}}=\frac{{J}_{i}^{y}}{4}-\frac{{J}_{i}^{z}{J}_{i}^{y}{J}_{i}^{z}}{3},$$

(22)

$${\tilde{S}}_{i,xy}^{z}={P}_{\frac{3}{2}}\,{S}_{i,xy}^{z}\,{P}_{\frac{3}{2}}=\frac{3{J}_{i}^{z}}{4}-\frac{{J}_{i}^{z}{J}_{i}^{z}{J}_{i}^{z}}{3},$$

(23)

$${\tilde{n}}_{i,xy}={P}_{\frac{3}{2}}\,{n}_{i,xy}\,{P}_{\frac{3}{2}}=\frac{3}{4}-\frac{{({J}_{i}^{z})}^{2}}{3}.$$

(24)

Unlike the simple Heisenberg model that only involves the linear spin operators, the effective model involves the spin products with two or three “J” operators. These operators are high order magnetic multipolar moments and are able to switch the local spin state from one J state to any other states, and thus quantum fluctuations are strongly enhanced. In terms of the “J” operators, the effective model is rather difficult to be tackled with, and the conventional wisdom cannot provide more physical intuition. Instead, we turn to the perspective of the emergent Kugel-Khomskii physics where the previous knowledge and theoretical techniques about the Kugel-Khomskii model can be adapted^{3,4,12,108,109,110}. For this purpose, we merely need to show the Kugel-Khomskii structure and reduce the effective model into the standard Kugel-Khomskii form.

We first make the following mapping between the effective spin states and the fictitious spin and orbital states,

$$\left\vert {J}_{i}^{z}=+\frac{3}{2}\right\rangle \equiv \left\vert {s}_{i}^{z}=+\frac{1}{2},{\tau }_{i}^{z}=+\frac{1}{2}\right\rangle ,$$

(25)

$$\left\vert {J}_{i}^{z}=+\frac{1}{2}\right\rangle \equiv \left\vert {s}_{i}^{z}=+\frac{1}{2},{\tau }_{i}^{z}=-\frac{1}{2}\right\rangle ,$$

(26)

$$\left\vert {J}_{i}^{z}=-\frac{1}{2}\right\rangle \equiv \left\vert {s}_{i}^{z}=-\frac{1}{2},{\tau }_{i}^{z}=-\frac{1}{2}\right\rangle ,$$

(27)

$$\left\vert {J}_{i}^{z}=-\frac{3}{2}\right\rangle \equiv \left\vert {s}_{i}^{z}=-\frac{1}{2},{\tau }_{i}^{z}=+\frac{1}{2}\right\rangle .$$

(28)

To distinguish from the physical spin “S”, we use the little “s” to label the fictitious spin, and use “τ” to label the fictitious orbital. Although τ is labeled as the “orbital”, the transformation under the time reversal differs from the usual orbital degree of freedom. This is because ${{{\mathcal{T}}}}\left\vert {J}_{i}^{z}=m\right\rangle ={i}^{2m}\left\vert {J}_{i}^{z}=-m\right\rangle$ and ${{{\mathcal{T}}}}$ does not modify the orbital degree of freedom for the usual real orbital wavefunctions. Here, the $\left\vert {J}_{i}^{z}=3/2\right\rangle$ and $\left\vert {J}_{i}^{z}=1/2\right\rangle$ have a different sign structure under the time reversal operation. This is the case even if we switch the assignment in Eq. (27) and Eq. (28). With the assignment in the above equations, the original J related operators can be expressed as

$${\tilde{S}}_{i,xy}^{x}=\frac{1}{3}{s}_{i}^{x}\,(1-2{\tau }_{i}^{z}),$$

(29)

$${\tilde{S}}_{i,xy}^{y}=\frac{1}{3}{s}_{i}^{y}\,(1-2{\tau }_{i}^{z}),$$

(30)

$${\tilde{S}}_{i,xy}^{z}=\frac{1}{3}{s}_{i}^{z}\,(1-2{\tau }_{i}^{z}),$$

(31)

$${\tilde{n}}_{i,xy}=\frac{1}{3}(1-2{\tau }_{i}^{z}).$$

(32)

Collecting all the interactions, we obtain the emergent Kugel-Khomskii model, and the interaction on the other bonds can be generated likewise. This interaction now carries the basic features of the conventional Kugel-Khomskii model with the following expression,

$${H}_{{{{\rm{eff}}}}}^{z}=\mathop{\sum}\limits_{\langle ij\rangle }\frac{4J}{9}\left({{{{\boldsymbol{s}}}}}_{i}\cdot {{{{\boldsymbol{s}}}}}_{j}+\frac{1}{4}\right)\left(\frac{1}{2}-{\tau }_{i}^{z}\right)\left(\frac{1}{2}-{\tau }_{j}^{z}\right).$$

(33)

The couplings on the remaining bonds can be obtained by the symmetry transformation. The understanding from the Kugel-Khomskii physics can then be applied. As the model from the J = 3/2 Mott insulators can be regarded as the parent model for many anisotropic spin-1/2 models, one can show that, with the strong easy-plane single-ion anisotropy like $+D{({J}_{i}^{x}+{J}_{i}^{y}+{J}_{i}^{z})}^{2}$, the effective spin-1/2 model carries a Kitaev-like interaction in addition to other anisotropic interactions.

Remarks for and beyond Kugel-Khomskii physics

We here give a few remarks based on the emergent Kugel-Khomskii physics for the multiflavor Mott insulators in quantum materials. The identification and separation of the effective spin and orbital degrees of freedom in these systems already indicates the distinct physical properties of these effective degrees of freedom. For instance, the effective spin is odd under time reversal and is usually related to the magnetic moment, while the effective orbital is often even under time reversal. In the case of the ordering phenomena, the effective spin and the effective orbital do not have to order at the same time^111,112,113, which is already feasible from the Ginzburg-Landau sense⁴⁰. Thus, one could have separate ordering transitions with different transition temperatures for the effective spin and orbital, respectively.

Due to the discrete lattice symmetry, the orbital ordering is often associated with the change of the crystal structure, leading to the lattice distortion and the phonon anomaly^112,113. This can be a useful experimental diagnosis of the effective orbital ordering. Like the traditional Kugel-Khomskii physics, the orbital order would reconstruct the spin model, the orientations of the orbital ordering then give rise to the spatially-anisotropic spin models that are often quite rare for the conventional magnets^114,115,116. For example, the quasi-one-dimensional Heisenberg chain physics in KCuF₃ was expected to arise from the orbital order that enhances the spin coupling along one lattice direction¹¹⁴. Such dimensional reduction is particularly useful to link the low-dimensional physics with the high-dimensional models and systems. The one-dimensional spin-1 Haldane-chain physics (and the symmetry protected topological order) was invoked for FeSe upon nematic or orbital ordering⁶⁹ and certain spin-1 honeycomb lattice Kugel-Khomskii model¹¹⁷.

Even in the ordered phases, the elementary excitations for these multiflavor Mott insulators are fundamentally different from the conventional magnets in the large-S limit. The latter is the simple spin-wave excitation. In contrast, the former is characterized by many branches of spin-wave-like excitations. For the conventional magnets with the large-S moment, the simple quadratic spin interaction merely creates the magnon excitation with the ±1 change of the spin quantum number and thus supports one spin-wave branch for each magnetic sublattice. In the multiflavor Mott insulators, the effective model could excite the system from the ordered flavor of state to all other flavor of states if the system is simply ordered to one flavor of state⁴⁰. For the emergent Kugel-Khomskii model with the spin-1/2 and the pseudospin-1/2 moments, there would be three branches of spin-wave-like excitations for each sublattice⁶⁰. These excitations are sometimes referred as the flavor-wave excitations.

It was previously argued that, the more exotic quantum states could be stabilized by the enhanced quantum fluctuations due to the enlarged local Hilbert space and the involved multiflavor interactions for the multiflavor Mott insulators. It has been shown that, the flavor ordered stripe state¹¹⁸, the valence bond crystal^119,120,121, Dirac spin liquid^45,122, spinon Fermi surface state¹²³ and chiral spin liquids^{56,119,120,122,124} could be stabilized by the enhanced quantum fluctuations and the large symmetry group like SU(N). The candidate systems are Ba₃CuSb₂O₉^125,126,127, the J = 3/2 Mott insulators like ZrCl₃^45,128, the double moiré layers with transition metal dichalcogenides or graphenes, as well as the ultracold atom contexts with the more exact symmetries that will be discussed later. The double moiré layers with transition metal dichalcogenides or graphenes have the enlarged symmetry like SU(4) or SU(8) from their spin, valley and layer freedom^124,129,130. Some of the results for the multiflavor Mott insulator of this context overlap with the ones from the ultracold atoms, and we do not give much discussion here.

In general, the multiflavor Mott insulator provides a feasible platform to study the exotic quantum states. To manipulate and probe these exotic quantum states and excitations, one necessarily has to understand the role of different degrees of freedom and how they are coupled to the external probes such as electric field, magnetic field, strains, electric currents and so on. The identification and separation of the degrees of freedom becomes much more meaningful for this practical purpose. For example, the selective coupling with the spin/magnetic excitations instead of the orbital modes or non-magnetic excitations in the neutron scattering and NMR measurements is crucial to understand the excitation structures for the conventional Kugel-Khomskii systems and other multiflavor Mott insulators¹³¹. The presence of the layer degrees of freedom in the multiflavor Mott insulators for the double moiré layers with transition metal dichalcogenides or graphenes allows the design of dragging type of transport measurements for the interlayer transport behaviors in the exotic quantum phases¹³⁰.

SU(N) and Sp(N) magnetism with ultracold fermions

Fermionic systems with large internal degrees of freedom not only exist in quantum material as multiflavor Mott insulators, but also include large-spin alkali and alkaline-earth fermions. These atoms often exhibit large hyperfine spins, whose 2S + 1 hyperfine atomic levels certainly form a high representation of the SU(2) group. As we have mentioned in the Introduction, theoretically it was proposed to study such systems from the perspective of exact high symmetries such as SU(N) and Sp(N)^47,49,132, in which the high symmetries greatly enhance the connectivity among different internal components. Consequently, they share very similar features in quantum magnetism in their Mott-insulating states to the Kugel-Khomskii systems despite the huge difference between their energy scales^{47,48,49,50,51,52,55,56,124,133,134,135,136}. The unprecedented exact high symmetries of the ultracold alkaline and alkaline-earth atom systems allow us to access many unconventional theoretical models and quantum phases that are a bit difficult to realize in quantum materials and only occur under the theoretical idealization. Unlike the previous discussion about the multiflavor Mott insulators and the Kugel-Khomskii physics in quantum materials, our review below will explore the models and the consequences from the high symmetries in the Mott insulating regime for the ultracold atoms.

For the alkaline-earth atoms, their hyperfine-spins completely come from the nuclear spins since their atomic shells are filled. Interactions from atomic scatterings are independent from spin components since the nuclei are deeply inside. This is the microscopic reason of their SU(N) interactions. Nevertheless, for the alkali high spin atoms, generally speaking, their interactions are spin-dependent. In this case, the SU(N) symmetry is not generic, and the Sp(N) symmetry is the next highest^{47,48,51,52,55}. As a concrete example, below we will use the simplest case of large spin ultra-cold fermions, i.e., spin-3/2 fermions (e.g. ¹³²Cs, ⁹Be, ¹³⁵Ba, ¹³⁷Ba, ²⁰¹Hg). It can be proved that such systems possess an exact and generic symmetry of Sp(4), or, isomorphically SO(5). Various many-body orders including quantum magnetism in different tensor channels, pairing and density-wave orders^47,48. The generic Sp(4) symmetry here plays the role of the SU(2) symmetry in the spin-1/2 systems.

Such a high symmetry group of the ultracold atom systems requires a fine-tuning of the parameters for the Kugel-Khomskii model in quantum materials. For instance, the symmetry can be tuned to SU(4) for the e_g orbitals. In this limit, the Kugel-Khomskii model in the quantum materials is connected to the SU(N = 4) model in the ultracold atom context in the narrow sense of the model relation. For the t_2g orbitals in quantum materials, more complicated operators are involved in the Kugel-Khomskii model⁷, and less symmetries are explicitly present. On the contrary, there were efforts⁴⁹ that introduce two atomic orbitals for the alkaline earth atoms such that the resulting model is the SU(N > 2) extension of the original the Kugel-Khomskii model where the SU(2) spin is replaced by the SU(N > 2) spin while the orbital sector of the electrons is replaced by the atomic orbitals on the optical lattices.

The “baryon-like” quantum magnetism is further reviewed. In quantum chromodynamics (QCD), quarks form the fundamental representation of the SU(3) group. Three quarks of all the R, G, B components form a baryon, a color singlet, while, two quarks cannot form a color singlet. Similarly, the SU(N) fermion systems are characterized by the N-fermion correlations, which leads to the multi-fermion SU(N) singlet clustering ordering^134,137,138. Similarly, the physics of an SU(N) quantum magnet can also be beyond the SU(2) case which is often studied condensed matter systems. In spite of the huge difference of energy scales, the large-spin cold fermions could exhibit similar physics to that in QCD^{133,135,139,140}.

The spin-3/2 Hubbard model–the generic Sp(4) symmetry

We generalize the spin-1/2 Hubbard model to spin-3/2 case. Only the standard spin SU(2) symmetry is assumed at the beginning, and the hidden Sp(4) symmetry will be explained. The free-fermion part H₀ reads,

$${H}_{0}=-t\mathop{\sum}\limits_{\langle ij\rangle ,\sigma }\left({\psi }_{i\sigma }^{{\dagger} }{\psi }_{j\sigma }+h.c.\right)-\mu \mathop{\sum}\limits_{i}n(i),$$

(34)

where ψ_σ is the 4-component fermion spinor operator and σ represents the eigenvalues of s_z with σ = ± 3/2, ± 1/2; 〈ij〉 denotes the nearest neighboring bonding; $n={\psi }_{\sigma }^{{\dagger} }{\psi }_{\sigma }$ is the particle-number operator. If two spin-3/2 fermions are put on the same site, their total spin is either 0 (singlet) or 2 (quintet). Then the onsite interactions are represented in the spin SU(2) language as

$${H}_{int}={U}_{0}\mathop{\sum}\limits_{i}{P}_{0}^{{\dagger} }(i){P}_{0}^{}(i)+{U}_{2}\mathop{\sum}\limits_{i,-2\le m\le 2}{P}_{2m}^{{\dagger} }(i){P}_{2m}^{}(i),$$

(35)

where ${P}_{0}^{{\dagger} }$ and ${P}_{2,m}^{{\dagger} }$ are the pairing operators in the singlet and quintet channels, respectively. Their interaction strengths are denoted as U₀ and U₂, respectively.

In fact, the Hamiltonian of Eqs. (34) plus (35) possesses an Sp(4) symmetry, which is the double covering group of SO(5). For simplicity, we will use Sp(4) and SO(5) interchangeably below. The Sp(4) algebra is constructed below. In addition to spin and charge, spin quadrupole and octupole operators are necessary in the spin-3/2 system. As shown in Methods, the five spin quadrupole matrices are actually the Dirac Γ-matrices denoted as Γ^a(1 ≤ a ≤ 5), and their commutators are defined as ${\Gamma }^{ab}=\frac{i}{2}[{\Gamma }^{a},{\Gamma }^{b}](1\le a \,<\, b\le 5)$. The spin-quadrupole operators n_a are defined as

$${n}_{a}=\frac{1}{2}{\psi }_{\alpha }^{{\dagger} }{\Gamma }_{\alpha \beta }^{a}{\psi }_{\beta },\quad (1\le a\le 5).$$

(36)

On the other hand, spin and spin octupole operators, in total 10 operators, can be reorganized as^47,132

$${L}_{ab}=-\frac{1}{2}{\psi }_{\alpha }^{{\dagger} }{\Gamma }_{\alpha \beta }^{ab}{\psi }_{\beta }.$$

(37)

L_ab’s span the Sp(4), or, isomorphically, SO(5) group, and n_a’s transform as a 5-vector of SO(5). L_ab and n_a together span the algebra of the SU(4) group, which is isomorphic to SO(6) group. Hence, Sp(4) is a subgroup of SU(4). Among them, three diagonal operators commute with each other:

$$\begin{array}{rcl}{L}_{15}=\frac{1}{2}({n}_{\frac{3}{2}}+{n}_{\frac{1}{2}}-{n}_{-\frac{1}{2}}-{n}_{-\frac{3}{2}}),\\ {L}_{23}=\frac{1}{2}({n}_{\frac{3}{2}}-{n}_{\frac{1}{2}}+{n}_{-\frac{1}{2}}-{n}_{-\frac{3}{2}}),\\ {n}_{4}=\frac{1}{2}({n}_{\frac{3}{2}}-{n}_{\frac{1}{2}}-{n}_{-\frac{1}{2}}+{n}_{-\frac{3}{2}}).\end{array}$$

(38)

Now it is ready to reveal the hidden Sp(4) symmetry. ψ_α is a 4-component spinor of SU(4) and n is an SU(4) scalar. Hence, H₀ is in fact SU(4) invariant, of course, Sp(4) invariant. It is enlightening to formulate Eq. (35) as

$${H}_{int}=\frac{3{U}_{0}+5{U}_{2}}{16}\mathop{\sum}\limits_{i}{\left[n(i)-2\right]}^{2}+\frac{{U}_{0}-{U}_{2}}{4}\mathop{\sum}\limits_{i,1\le a\le 5}{n}_{a}^{2}(i).$$

(39)

Since the n^a operators form a 5-vector, Eq. (39) is also Sp(4) invariant.

The superexchanges at quarter-filling

Under sufficiently strong repulsive interactions, spin-3/2 Hubbard model will enter the Mott-insulating state even at 1/4-filling, i.e., one fermion per site. In principle, the magnetic superexchange interaction can be expressed in terms of the bi-linear, bi-quadratic and bi-cubic SU(2) Heisenberg terms of S_i ⋅ S_j, ${({{{{\bf{S}}}}}_{i}\cdot {{{{\bf{S}}}}}_{j})}^{2}$, ${({{{{\bf{S}}}}}_{i}\cdot {{{{\bf{S}}}}}_{j})}^{3}$, respectively. Nevertheless, the superexchange exists in the bond-spin singlet and quintet channels, respectively. Hence, only two of these three terms are independent, which corresponds to the Sp(4) symmetric case.

The above situation is similar to the spin-1 Heisenberg model with both bi-linear and bi-quadratic terms. Its general form is controlled by the parameter angle θ,

$${H}_{spin1}=J\mathop{\sum}\limits_{\langle ij\rangle }\left(\cos \theta \,({{{{\bf{S}}}}}_{i}\cdot {{{{\bf{S}}}}}_{j})+\sin \theta \,{({{{{\bf{S}}}}}_{i}\cdot {{{{\bf{S}}}}}_{j})}^{2}\right).$$

(40)

Eq. (40) exhibits two different types of SU(3) symmetries, a “uniform” one and a “staggered” one. The “uniform” one appears at θ = π/4 and 5π/4, i.e., every site lies in the fundamental representation. The “staggered” one appears at θ = π/2 and 3π/2, in which two sublattices lie in the fundamental and anti-fundamental representations, respectively. They are explained in Methods.

It is enlightening to express the superexchange model of spin-3/2 fermions in the explicitly Sp(4) invariant form¹³³,

$${H}_{ex}=\mathop{\sum}\limits_{\langle ij\rangle }\left\{{J}_{L}\mathop{\sum}\limits_{1\le a < b\le 5}{L}_{ab}(i){L}_{ab}(j)+{J}_{N}\mathop{\sum }\limits_{a=1}^{5}{n}_{a}(i){n}_{a}(j)\right\}.$$

(41)

Each site lie in the fundamental representation of Sp(4), i.e., a 4-component spinor. J_L = (J₀ + J₂)/4 and J_N = (3J₂ − J₀)/4, and J₀ and J₂ are exchange strength in the singlet and quintet channels, respectively. At the level of the 2nd order perturbation theory, they are

$${J}_{0}=4{t}^{2}/{U}_{0},\quad {J}_{2}=4{t}^{2}/{U}_{2}.$$

(42)

Sp(4) is a rank-2 Lie algebra. Hence, it has three good quantum numbers, i.e., one more compared to those of SU(2). They are defined as

$$C=\mathop{\sum}\limits_{1\le a < b\le 5}{\left\{\mathop{\sum}\limits_{i}{L}_{ab}(i)\right\}}^{2},$$

(43)

$${L}_{15}^{tot}=\mathop{\sum}\limits_{i}{L}_{15}(i),\quad {L}_{23}^{tot}=\mathop{\sum}\limits_{i}{L}_{23}(i).$$

(44)

C is the Sp(4) Casimir in analogy to total spin square of an SU(2) system; ${L}_{15}^{tot}$ and ${L}_{23}^{tot}$ are the analogies to the total S_z.

The uniform and staggered SU(4) symmetries

Similar to the spin-1 model of Eq. (40), Eq. (41) exhibits two different types of SU(4) symmetries. The first case takes place at J₀ = J₂, or, equivalently, U₀ = U₂, denoted as SU(4)_A below. Eq. (41) is reduced to the SU(4) Heisenberg model

$${H}_{A}=J\mathop{\sum}\limits_{\langle i,j\rangle }\left\{{L}_{ab}(i){L}_{ab}(j)+{n}_{a}(i){n}_{a}(j)\right\},$$

(45)

for which each site lies in the fundamental representation of SU(4), and J = J₀/2 = J₂/2. H_A is equivalent to the SU(4) Kugel-Khomskii type model^{7,139,141,142,143} as reviewed in the previous parts, which is similar to the “uniform” case of SU(3).

The second SU(4) symmetry is similar to the “staggered” case of SU(3), which takes place in a bipartite lattice and in the limit of U₂ → + ∞, i.e., J₂ = 0. To see this explicitly, the particle-hole transformation ${\psi }_{\alpha }\to {R}_{\alpha \beta }{\psi }_{\beta }^{{\dagger} }$ is performed to one sublattice, and the other sublattice is left unchanged, where R is the charge conjugation matrix

$$R=\left(\begin{array}{cc}0&i{\sigma }_{2}\\ i{\sigma }_{2}&0\\ \end{array}\right).$$

(46)

Under this operation, the fundamental representation of SU(4) transforms to its anti-fundamental representation whose Sp(4) generators and vectors become ${L}_{ab}^{{\prime} }={L}_{ab}$ and ${n}_{a}^{{\prime} }=-{n}_{a}$. The Sp(4) generators remain invariant under the particle-hole transformation is due to its pseudo-reality. Then Eq. (41) can be recast to

$${H}_{B}={J}^{{\prime} }\mathop{\sum}\limits_{\langle i,j\rangle }\left\{{L}_{ab}^{{\prime} }(i){L}_{ab}^{}(j)+{n}_{a}^{{\prime} }(i){n}_{a}(j)\right\},$$

(47)

where ${J}^{{\prime} }={J}_{0}/4$. Eq. (47) is SU(4) invariant again.

These two SU(4) symmetries imply fundamentally different physics. Both of them have high energy analogies: The “uniform” SU(4)_A is baryon-like, while that of SU(4)_B is meson-like. In the former case, every site belongs to the fundamental representation. At least 4 sites are required to form an SU(4) singlet with the wavefunction represented by

$$\frac{1}{\sqrt{4!}}{\epsilon }_{\alpha \beta \gamma \delta }{\psi }_{\alpha }^{{\dagger} }(1){\psi }_{\beta }^{{\dagger} }(2){\psi }_{\gamma }^{{\dagger} }(3){\psi }_{\delta }^{{\dagger} }(4)\left\vert \Omega \right\rangle ,$$

(48)

where 1 ~ 4 are site indices, ϵ_αβγδ is the rank-4 fully antisymmetric tensor, and Ω is the vacuum. Hence, dramatically different from the SU(2) case, two sites across a bond cannot form a singlet. Quantum magnetism based on H_A exhibits strong features of the 4-site correlation, i.e., a baryon-like state.

For the case of SU(4)_B, two sites are sufficient to form an SU(4) singlet via the charge conjugation matrix as

$$\frac{1}{2}{R}_{\alpha \beta }{\psi }_{\alpha }^{{\dagger} }(1){\psi }_{\beta }^{{\dagger} }(2)\left\vert \Omega \right\rangle .$$

(49)

This can be viewed as a large-N version of the usual spin-1/2 SU(2) Heisenberg model.

Quantum magnetism of an Sp(4) spin chain

The Sp(4) quantum magnetism occurs in the strong repulsive interaction regime at 1/4-filling, i.e., one particle per site. For the 1D system, the bosonizaiton study shows that the charge gap opens due to the relevancy of the 4k_f-Umklapp process at the Luttinger parameter 0 < K_c < 1/2^132,134. Then the low energy physics is captured by the Sp(4) superexchange process described by Eq. (41).

The 1D phase diagram is sketched in Fig. 5. A quantum phase transition takes place at the SU(4) symmetric point of J₀ = J₂, or, U₀ = U₂: When J₀ > J₂ (U₀ < U₂), the system develops a spin gap with the presence of the spin Peierls distortion, while at J₀ ≤ J₂ (U₀ ≤ U₂), the system enters a gapless spin liquid phase and maintains the translation symmetry¹³². The nature of this transition is Kosterlitz-Thouless like. For convenience, a parameter angle θ is employed to represent J_0,2 as

$${J}_{0}=\sqrt{2}\cos \theta ,{J}_{2}=\sqrt{2}\sin \theta .$$

(50)

Hence, the gapless spin liquid phase lies at 0 ≤ θ ≤ 45^∘, while the spin Peirerls phase exhibiting dimerization lies at 45^∘ < θ ≤ 90^∘.

Fig. 5: Phase diagram of the 1D Sp(4) spin chain in terms of J₀ and J₂ with $\tan \theta ={J}_{0}/{J}_{2}$.

The SU(4)_A and SU(4)_B symmetries are located along the lines of θ = 45^∘ and 90^∘, respectively. The dimerized spin Peierls phase appears at 90^∘ ≥ θ > 45^∘, and the gapless spin liquid phase is located at 45^∘ ≥ θ ≥ 0^∘. Figure is adapted from ref. ¹³⁶.

The above field-theoretical analysis was confirmed by the density-matrix-renormalization-group (DMRG) simulations¹³⁶. The ground state spin gap Δ_sp is defined as the energy difference between the ground state and the lowest Sp(4) multiplet. The DMRG results show that for the cases of θ > 45^∘, i.e., J₂/J₀ < 1, Δ_sp’s saturate as enlarging the lattice size, which signatures the opening of spin gap. On the other hand, Δ_sp’s vanish at θ ≤ 45^∘ demonstrating gapless ground states. These results show that the phase boundary is located at θ = 45^∘ with the SU(4)_A symmetry, on which the system is also gapless.

The spin gapped and gapless phases exhibit different characteristic oscillations. The DMRG results of the nearest neighbor (NN) correlations of the Sp(4) generators are calculated with the open boundary condition. As an example, 〈L₁₅(i)L₁₅(i + 1)〉 is shown in Fig. 6 and those of other Sp(4) generators are the same. In the spin gapped phase, say, at θ = 60^∘, a characteristic 2-site oscillation, corresponding to a dimer pattern, is pinned by the open boundary condition. It does not decay as moving to the center of the chain implying the presence of long-range-ordering in the thermodynamic limit. This is consistent with the 4k_f charge-density-wave ordering, exhibiting as dimerization. In contrast, in the gapless region 0 ≤ θ ≤ 45^∘, the open boundary induces a power-law decay. For example, at θ^∘ = 30^∘, it approximately exhibits a 4-site oscillation, and the periodicity is the same for 0 ≤ θ ≤ 45^∘. This is in agreement with the dominant 2k_f spin correlations in the bosonization analysis.

Fig. 6: The DMRG results of the nearest neighbor correlations of 〈L₁₅(i)L₁₅(i + 1)〉 under the open boundary condition.

The parameter values are θ = 60^∘ in (a) and θ = 30^∘ in (b). A 2-site periodicity appears in (a) and a 4-site periodicity with a power-law decay shows up in (b). Figures are adapted from ref. ¹³⁶.

In the gapless phase, i.e., 0 ≤ θ≤ 45^∘, the decay powers of the two-point correlation functions are discussed below. The correlation can be expressed asymptotically as

$$\langle O({i}_{0})O(i)\rangle \propto \frac{\cos \frac{\pi }{2}| {i}_{0}-i| }{| {i}_{0}-i{| }^{\kappa }},$$

(51)

where the operator O = L, or, n represents an operator in the 10-generator channel, or, the 5-vector channel of the Sp(4) group, respectively. The critical exponents for L and n along the SU(4)_A line (θ = 45^∘) should be the same, as fitted by κ ≈ 1.52, which is in a good agreement with the value of 1.5 from the bosonization analysis^134,144,145. As away from the SU(4)_A line, the degeneracy between L and n is lifted. The simulations show that the critical exponent for L shows κ_L < 1.5, while that for n exhibits κ_n > 1.5.

The 2D Sp(4) magnetism – Exact diagonalization results

The 2D Sp(4) magnetism based on Eq. (41) is very challenging, whose global phase diagram remains unsolved. Especially, in the region of 0 < θ ≤ 45^∘, the phase is dominated by the baryon-type physics which is quite different from the SU(2) quantum magnetism^133,142. On the other hand, the physics along the SU(4)_B line, i.e, θ = 90^∘, is relatively clear. Quantum Monte Carlo simulations (QMC) show a long-range Néel ordering with a much weaker Néel moment compared to the SU(2) case¹⁴⁶. However, except this special case, the sign-problem appears and no conclusive results are known. Along the SU(4)_A line, it is likely that there exist a strong plaquette-like correlation as proposed by Bossche et. al.^139,140. Such a state is a baryon-like state, exhibiting the plaquette type correlations, and each plaquette is a 4-site SU(4) singlet state.

Below we review the ground state correlations via the exact diagonalization (ED) method. Even though it can only be applied to a small 4 × 4 cluster, the associated ground state profiles at different values of θ still yield valuable information to speculate the thermodynamic limit. Based on the ED result, a phase transition would be expected as follows: A Néel order-like state breaking the Sp(4) symmetry appears at large values of θ, in particular close to 90^∘; while an Sp(4) singlet ground state breaking the translation symmetry exits as lowering θ to smaller values.

We use structure factors to describe the ground state properties. A structure factor converges to a finite value in the thermodynamic limit in the presence of long-range ordering^146,147. Nevertheless, for a small cluster, the tendency of ordering is reflected as the peak of a structure factor at a characteristic momentum. The following structure factors will be employed,

$${S}_{O}({{{\bf{q}}}})=\frac{1}{N}\mathop{\sum}\limits_{i,j}{e}^{i{{{\bf{q}}}}\cdot ({{{{\bf{r}}}}}_{i}-{{{{\bf{r}}}}}_{j})}\langle G| O(i)O(j)| G\rangle ,$$

(52)

where O = L referring to an operator of the Sp(4) generators, and O = n referring to an operator of the Sp(4) vectors. For later convenience, we use the following symbols to represent the crystal momenta Γ = (0, 0), X = (π, 0), and M = (π, π).

Previous QMC simulations show the Néel long-range order of the SU(4) Heisenberg model defined in the fundamental and anti-fundamental representations in two sublattices¹⁴⁶. It corresponds to the SU_B(4) case with θ = 90^∘, which ensures the relation of

$${S}_{n}({{{\bf{q}}}})={S}_{L}({{{\bf{q}}}}+M),$$

(53)

due to the staggered definition of Sp(4) vectors n_a in Eq. (47). Back to the spin language, as θ = 90^∘ the classic energy can be minimized by arranging the two components of S_z = ± 3/2, or, of the other two components of S_z = ± 1/2 in a staggered way. These different configurations are equivalent under the Sp(4) transformations.

The SU_B(4) relation still approximately holds when θ is close to 90^∘. The Néel type correlation of L_ab also extends to a finite regime as θ deviates away from 90^∘, and in the same regime, n_a exhibits the uniform correlations. When 90^∘ ≥ θ ≳ 60^∘, the structure factor S_L(q) of the 10-generator channel peaks at the M-point, i.e., the nesting wavevector (π, π), as shown in Fig. 7(a). On the other hand, S_n(q) of the 5-vector channel is shown in Fig. 7(b), which peaks at the Γ-point, roughly in the same parameter range that S_L(q) develops a peak at the M-point.

Fig. 7: The structure factors for the 4 × 4 cluster with (a) S_L(q) in the Sp(4) generator sector and (b) S_n(q) in the Sp(4) vector sector.

The expressions of S_L and S_n are defined in Eq. (52). The inset in (a) is the comparison between S_L(π, 0) and S_L(π, π) versus θ. Figures are adapted from ref. ¹³⁶.

In contrast, at θ ≲ 60^∘, S_L(q) distributes relatively smooth over all momenta without dominant peaks. The distribution of S_n(q) also becomes smooth in a similar way to that of S_L(q). Nevertheless, at small values of θ < 18^∘, it becomes to peak at the M-point, i.e., at the momentum (π, π). In the case of SU(4)_A, i.e., θ = 45^∘, the 5-vector and the 10-generator channels become degenerate, hence, S_n(q) = S_L(q) for each q. As θ deviates away from 45^∘, S_L and S_n evolve differently.

To test the possibility towards a dimerized ground state, the susceptibilities of translation and rotational symmetry breakings are defined below. Two perturbations are added to the Hamiltonian Eq. (41), i.e.,

$${\hat{O}}_{dim}({{{\bf{Q}}}})=\mathop{\sum}\limits_{i}\cos ({{{\bf{Q}}}}\cdot {{{{\bf{r}}}}}_{i}){H}_{ex}(i,i+\hat{x}),$$

(54)

$${\hat{O}}_{rot}=\mathop{\sum}\limits_{i}[{H}_{ex}(i,i+\hat{x})-{H}_{ex}(i,i+\hat{y})].$$

(55)

The former breaks the translation symmetry and the latter breaks the rotation symmetry.

The ED results of susceptibilities associated with ${\hat{O}}_{dim}$ and ${\hat{O}}_{rot}$ show that both of them exhibit a peak in the interval of 60^∘ < θ < 70^∘136. Although no real divergences exist due to the finite size, sharp peaks would imply the tendency of long-range ordering. Hence, the results imply a tendency to break both translational and rotational symmetries, which is consistent with a columnar dimerization in the thermodynamic limit. In contrast, the plaquette ordering does not break the 4-fold rotational symmetry, which is not favored in this regime.

The tendency of the plaquette type ordering, i.e., the SU(4) analogy to the baryon state, shows up at small values of θ, i.e., θ < 50^∘ ~ 60^∘. To characterize such a state, the plaquette Sp(4) Casimir centered at r is defined as

$$C({{{\bf{r}}}})=\langle G| \mathop{\sum}\limits_{1\le a < b\le 5}{\left\{\mathop{\sum}\limits_{i}{L}_{ab}(i)\right\}}^{2}| G\rangle ,$$

(56)

where i runs over the four sites of this plaquette. A similar method was used for the SU(2) case to characterize the competing dimer and plaquette orders¹⁴⁸.

The relation of C versus θ is calculated under the open boundary condition. Based on the symmetry analysis, the values of C for three non-equivalent plaquettes A, B, and C are shown in Fig. 8. C(r_A) of the corner plaquette is significantly about one order smaller than C(r_B) and C(r_C) of the other two plaquettes at small values of θ. This suggests that an Sp(4) plaquette tendency is pinned down by the open boundary. As θ increases above 60^∘, the contrast decreases which means that the plaquette-type pattern is weakened. It is likely that there exists a strong plaquette-like correlation at small values of θ but significantly beyond the SU(4)_A line. It covers the entire region of θ < 45^∘ and also extends slightly above 45^∘. Nevertheless, further larger scale calculations are necessary for a conclusion.

Fig. 8: The local Casimir C(r) defined in Eq. (56) versus θ.

The plaquettes (A, B, C) are depicted in the right.

Quantum plaquette model for the 3D SU(4) magnetism

The SU(4)_A plaquette state discussed above has been shown as the exact solution to the ground state of a 2-leg ladder model of spin-3/2 fermions¹³³. However, due to the geometric constraint, it cannot resonate and is a valence-bond-solid type state.

The resonating quantum plaquette model (QPM) model was constructed in three dimensions^135,149, which is analogous to the quantum dimer model for the SU(2) magnet^150,151. The quantum dimer model has a gauge theory description to the Rohksar-Kivelson (RK) Hamiltonian, which is a compact U(1) gauge theory¹⁵¹. QPM also has a similar description as reviewed below, which is a high order gauge field theory. Recently it has received considerable attention in the context of fracton physics¹⁵².

Consider a three-dimensional (3D) cubic lattice SU(4) model in the limit that each plaquette has a strong tendency to form a local SU(4) singlet. The effective Hilbert space is spanned by all the plaquette configurations. They are subject to the constraint that every site belongs to one and only one plaquette.

Each unit cube possesses three flippable configurations: the pairs of plaquettes of left and right, top and bottom, and front and back denoted as A, B and C in Fig. 9, respectively. A Rokhsar-Kivelson (RK) type Hamiltonian is constructed as¹⁵⁰:

$$\begin{array}{l}H\,=\,-t\mathop{\sum}\limits_{{{{\rm{each}}}}\,{{{\rm{cube}}}}}\left\{\left\vert A\right\rangle \left\langle B\right\vert +\left\vert B\right\rangle \left\langle C\right\vert +\left\vert C\right\rangle \left\langle A\right\vert +h.c.\right\}\\ \,\,\,\,\quad\quad+\,V\mathop{\sum}\limits_{{{{\rm{each}}}}\,{{{\rm{cube}}}}}\left\{\left\vert A\right\rangle \left\langle A\right\vert +\left\vert B\right\rangle \left\langle B\right\vert +\left\vert C\right\rangle \left\langle C\right\vert \right\},\end{array}$$

(57)

where t is assumed to be positive, and v/t is arbitrary. By defining

$$\left\vert {Q}_{1,2}\right\rangle =\left\vert A\right\rangle +{e}^{\pm i\frac{2}{3}\pi }\left\vert B\right\rangle +{e}^{\mp i\frac{2}{3}\pi }\left\vert C\right\rangle ,$$

(58)

Eq. (57) is reformulated as

$$\begin{array}{l}H\,=\,t\mathop{\sum}\limits_{{{{\rm{each}}}}\,{{{\rm{cube}}}}}\left(\left\vert {Q}_{1}\right\rangle \left\langle {Q}_{1}\right\vert +\left\vert {Q}_{2}\right\rangle \left\langle {Q}_{2}\right\vert \right)\\ \,\quad\quad+\,(V-2t)\mathop{\sum}\limits_{{{{\rm{each}}}}\,{{{\rm{cube}}}}}\left(\left\vert A\right\rangle \left\langle A\right\vert +\left\vert B\right\rangle \langle B+| C\rangle \left\langle C\right\vert \right).\end{array}$$

(59)

The RK point here corresponds to V = 2t.

Fig. 9: Three flippable configurations in one cube.

Figure is adapted from ref. ¹³⁵.

The ground state to Eq. (59) is the equal weight superposition of all the plaquette configurations within the same topological sector that can be connected by local flips. The classical Monte Carlo simulation shows that a crystalline order of resonating cubes at this RK point¹⁴⁹. At V/t < 2, the system favors flippable cubes. For example, as shown in Ref. [¹⁴⁹], at v/t ≪ − 1, the ground state exhibit the columnar ordering. At v > 2t, both terms in Eq. (59) are positive-definite. Hence, all the configurations without flippable cubes are the ground states. All the transitions between different phases are of the first order.

It is often useful to extract the low energy physics of strong correlated systems by mapping it to gauge theory models. In fact, the effective gauge theory of the QPM was constructed as a special one - a high order gauge theory.

Each square face is denoted by the location of its center: $i+\frac{1}{2}\hat{u}+\frac{1}{2}\hat{v}$, where $\hat{u}=\hat{x},\hat{y}$ and $\hat{z}$. Each face is associated with a number n and a strong local potential $U{({n}_{i+\frac{1}{2}\hat{\mu }+\frac{1}{2}\hat{\nu }}-\frac{1}{2})}^{2}$ in the limit of U → ∞, such that the low-energy sector has only either n = 1 corresponding to an SU(4) singlet occupation, or n = 0, otherwise. The constraint is that the summation of n over all the 12 faces connecting to the same site equals to 1. A rank-2 symmetric traceless tensor electric field is defined as

$${E}_{i,\mu \nu }=\eta (i)\left({n}_{i+\frac{1}{2}\hat{\mu }+\frac{1}{2}\hat{\nu }}-\frac{1}{2}\right),$$

(60)

where η(i) = ± 1 depends on the sublattice that i belongs to. The constraint is simply reduced to

$${\nabla }_{x}{\nabla }_{y}{E}_{xy}+{\nabla }_{y}{\nabla }_{z}{E}_{yz}+{\nabla }_{z}{\nabla }_{x}{E}_{zx}=5\eta (i),$$

(61)

where ∇ is the lattice derivative. The canonical conjugate variable of E_i,μν is the vector potential A_i,μν, satisfying

$$[{E}_{i,\mu \nu },{A}_{j,\rho \sigma }]=i{\delta }_{ij}({\delta }_{\mu \rho }{\delta }_{\nu \sigma }+{\delta }_{\mu \sigma }{\delta }_{\nu \rho }).$$

(62)

A_i,μν can be expressed by the phase variable ${\theta }_{i+\frac{1}{2}\hat{\mu }+\frac{1}{2}\hat{\nu }}$, which is the canonical conjugate variable to n as ${A}_{i,\mu \nu }=\eta (i)\,{\theta }_{i+\frac{1}{2}\hat{\mu }+\frac{1}{2}\hat{\nu }}$. Because E_μν only takes integer values, A_μν is an compact field with period of 2π.

The plaquette flipping process changes the plaquette occupations. For this purpose, ${e}^{i{A}_{j,\nu \sigma }}$ is employed which changes the eigenvalue of E_i,μν by 1 since

$$[{E}_{i,\mu \nu },{e}^{i{A}_{j,\nu \sigma }}]=({\delta }_{\mu \rho }{\delta }_{\nu \sigma }+{\delta }_{\mu \sigma }{\delta }_{\nu \rho }){e}^{i{A}_{j,\nu \sigma }}.$$

(63)

The flipping term is represented as

$$\begin{array}{l}{H}_{t}\,=\,-t\left[\cos ({\nabla }_{z}{A}_{xy}-{\nabla }_{x}{A}_{yz})+\cos ({\nabla }_{x}{A}_{yz}-{\nabla }_{y}{A}_{zx})\right.\\ \qquad\quad+\,\left.\cos ({\nabla }_{y}{A}_{zx}-{\nabla }_{z}{A}_{xy})\right].\end{array}$$

(64)

The associated gauge invariant transformation is

$${A}_{\mu \nu }\to {A}_{\mu \nu }+{\nabla }_{\mu }{\nabla }_{\nu }f,$$

(65)

where f an arbitrary scalar function. The low-energy Hamiltonian of the system can be written as

$$\begin{array}{l}H\,=\,{H}_{t}+\mathop{\sum}\limits_{{{{\rm{each}}}}\,{{{\rm{cube}}}}}\left\{U[{E}_{xy}^{2}+{E}_{yz}^{2}+{E}_{zx}^{2}]\right.\\ \qquad\quad+\,\left.V[{({\nabla }_{x}{E}_{yz})}^{2}+{({\nabla }_{y}{E}_{zx})}^{2}+{({\nabla }_{z}{E}_{zx})}^{2}]\right\},\end{array}$$

(66)

under the constraint of Eq. (61).

The above formulation is a compact high order gauge theory, and the non-local topological defect excitations is crucial to determine whether its ground state is gapless or gapped. The standard analysis is the mapping to a height model, such that topological defects are represented by the vertex operators¹³⁵. Due to the proliferation of topological defects, the system is generally gapped for the whole phase diagram of the quantum model, which corresponds to crystalline orders. It would be interesting to further explore the possiblity of the 3D power-law phases as the analog to the Kosterlitz–Thouless transition¹⁵².

Summary and future perspective

In summary, we have reviewed a class of Mott insulators, namely, the multiflavor Mott insulators, in both quantum materials and ultracold atom systems. In both cases, the local Hilbert states in each unit cell of the multiflavor Mott insulators consist of more than two flavors in contrast to the conventional spin-1/2 Mott insulators, and also differs fundamentally from the conventional magnets with large S moments. They bear similarities to the orbital-active Mott-insulators but typically do not demand the explicit orbital degeneracy. These include the breathing/cluster magnets, the J = 3/2 Mott insulators in various transition metal oxides and rare-earth systems, twisted moiré heterostructures, the ultracold atom fermion systems with large hyperfine-spins, and various other physical contexts.

The study of the multiflavor Mott insulators will certainly broaden the research scope of quantum magnetism, and further enrich the activity of exploring exotic quantum states of matter. As the consequence of the large local Hilbert space, this class of Mott insulators exhibit the common feature of enhanced quantum fluctuations. Instead of being viewed from the large-S perspective, the appropriate viewpoint should be large-N. The low-energy superexchange models typically go beyond the SU(2) large-S Heisenberg models, and bear similarities to the Kugel-Khomskii models. They are expected to exhibit a variety unusual physics including the spin-multipolar ordering, the “baryon-like” like physics, and even more exotic spin-liquid states. Quantum magnetism associated with large symmetries could exhibit certain similarities to QCD in spite of dramatically different energy scales. If each lattice site belongs to the fundamental representation of SU(N) with N > 2, the SU(N) singlet typically involves more than two sites leading to the “baryon-type” physics.

The review mainly focuses on the multiflavor Mott insulating states and the associated pairwise interactions with and without the enlarged symmetries in the strong Mott regime where the charge fluctuation is suppressed. In many cases, the system could be proximate to the Mott transition and in the weak Mott regime. In this regime, the charge fluctuation is strong, and the more appropriate description of the underlying physics would be in terms of the multiflavor Hubbard model or the exchange model involving the multiflavor ring exchange interactions. It is well-known that, the strong charge fluctuation of the weak Mott insulators or proximity to the Mott transition could stabilize various exotic quantum phases even for the spin-1/2 degrees of freedom^153,154. The multifavor Mott insulators actually bring a bit more interesting phenomena with their multiple favors of degrees of freedom. One well-known consequence is the Pomeranchuk effect or Pomeranchuk cooling that has been observed in magic-angle graphene^155,156 and been used to cool to much lower temperatures in the SU(N) alkaline-earth atoms^54,157. Although both the strong charge fluctuation and the multiflavor Hilbert space could drive to exotic quantum phases and phenomena, the exotic phases that are favored by them, however, can be different^123,158,159. There has not been much systematic study in this direction. Potentially this can be an interesting direction for both theories and experiments. In the ultracold atom systems, the strength of correlation can be effectively tuned experimentally to access both strong and weak Mott regimes. With the modern fabricating and controlling technology, one can effectively tune many physical properties of the transition metal dichalcogenide (moiré) heterostructures such as the correlation and doping, and drive the system through the Mott transition³³. In fact, the type of charge fluctuations in the early study of the Hubbard model is mostly concerned about the physical processes above the Mott gap that involve the double occupation¹⁵³. The existence of large numbers of correlated moiré systems and the partially-filled correlated materials supports the cluster localization for which the possibility of the sub-Mott-gap charge fluctuations can be relevant^20,22,29,160. Distinct charge fluctuations would bring rather distinct physics and consequences on the remaining spin and orbital degrees of freedom, especially since the charge sector often has a higher energy scale than the spin and orbital and thus may be considered first. The combination of distinct charge fluctuations and high symmetries is another direction to go. When the charge sector is incorporated as one extra flavor of the multiflavor Mott systems, not only the physics becomes richer from the interplay of more degrees of freedom, but also the scope of multiflavor Mott insulators can be further broadened. Thus, the multiflavor Mott systems and physics are not limited to the contents and the classification in the current review, and many more systems may be recast into this topic. We expect many exciting and new results to emerge in this field in the near future.

Methods

The SU(3) viewpoint to the spin-1 Heisenberg model with the bi-quadratic interaction

The two-band Hubbard model described by Eq. (1) exhibits different types of Mott-insulating states depending on the filling factor. The case of quarter-filling exhibiting the orbital degree of freedom is already discussed in the main text, whose low-energy physics is described by a variant of the SU(4) Kugel-Khomskii model. On the other hand, for the case of half-filling, i.e., two electrons per site, each orbital is occupied by one electron. Therefore, the orbital degree of freedom is quenched, and only the spin degree of freedom remains. The Hund’s rule coupling aligns two electron spins and forms the spin-1 moment.

In this case, the leading order super-exchange process across a bond is the switching of one pair of electrons, which is at the level of the second order perturbation theory. It gives rise to the standard bi-linear Heisenberg superexchange term as

$${H}_{bl}={J}_{2}\mathop{\sum}\limits_{\langle ij\rangle }{{{{\bf{S}}}}}_{i}\cdot {{{{\bf{S}}}}}_{j}.$$

(67)

Moreover, the fourth order perturbation theory yields the bi-quadratic superexchange term as

$${H}_{bq}={J}_{4}\mathop{\sum}\limits_{\langle ij\rangle }{\left({{{{\bf{S}}}}}_{i}\cdot {{{{\bf{S}}}}}_{j}\right)}^{2}.$$

(68)

Based on their perturbation orders, J₂ and J₄ are estimated at

$${J}_{2}\propto \frac{{t}^{2}}{U},\quad {J}_{4}\propto \frac{{t}^{4}}{{U}^{3}}.$$

(69)

Hence J₄ is much smaller than J₂. The other microscopic mechanism for the bi-quadratic interaction is the spin-lattice coupling. This happens when the direct exchange dominates and depends sensitively on the bond length. Integrating out the bonds, one can obtain the bi-quadratic spin interaction¹⁶¹. This has been applied to the Cr-based pyrochlores CdCr₂O₄ and HgCr₂O₄. As far as we are aware of, the sign of the bi-quadratic spin interaction is often negative and ferromagnetic-like.

The dominance of the bi-linear super-exchange process applies to even higher spin systems in quantum materials. The onsite SU(2) multiplets are denoted as $\left\vert S{S}_{z}\right\rangle$ with S ≥ S_z ≥ − S. The bi-linear superexchange can only connect these states via a one-dimensional weight diagram as

$$\left\vert S-S\right\rangle \leftrightarrow \left\vert S-S+1\right\rangle \ldots \leftrightarrow \left\vert SS\right\rangle .$$

(70)

Hence, spin flipping between $\left\vert S-S\right\rangle \leftrightarrow \left\vert SS\right\rangle$ takes 2S successive steps, which is at the 4S-th order perturbation theory. Therefore, no matter how large the total spin is, quantum spin fluctuations are basically a 1/S-effect.

The bi-quadratic term in Eq. (68) and bi-linear on in Eq. (67) can be treated at equal footing as controlled by the parameter angle θ as shown in Eq. (40). Although for physical systems of transition metal oxides, 0 < ∣θ∣ ≪ π/2, below the entire parameter space will be considered.

For one bond 〈ij〉, according to the SU(2) structure, the bond Hilbert space is divided into the singlet, triplet and quintet states. Their energies are denoted as E_s, E_t and E_q, respectively. It is straightforward to show that

$$\begin{array}{lll}{E}_{s}/J=-2\cos \theta +4\sin \theta ,\\ {E}_{t}/J=-\cos \theta +\sin \theta ,\\ {E}_{q}/J=\cos \theta +\sin \theta .\end{array}$$

(71)

At two special angles, degeneracy patterns appear beyond the SU(2) multiplet: A 6-fold degeneracy appears at θ = π/4 and 5π/4, and an 8-fold degeneracy appears at θ = π/2 and θ = 3π/2. They imply the hidden high symmetries of two different types of SU(3).

To show the SU(3) symmetry explicitly, the common representation of the Gell-mann matrix is adopted:

$$\begin{array}{l}{\lambda }_{1,ij}\,=\,{\delta }_{i1}{\delta }_{j2}+{\delta }_{i2}{\delta }_{j1},\quad {\lambda }_{2,ij}=-i{\delta }_{i1}{\delta }_{j2}+i{\delta }_{i2}{\delta }_{j1}\\ {\lambda }_{4,ij}\,=\,{\delta }_{i3}{\delta }_{j1}+{\delta }_{i1}{\delta }_{j3},\quad {\lambda }_{5,ij}=-i{\delta }_{i3}{\delta }_{j1}+i{\delta }_{i1}{\delta }_{j3}\\ {\lambda }_{6,ij}\,=\,{\delta }_{i2}{\delta }_{j3}+{\delta }_{i3}{\delta }_{j2},\quad {\lambda }_{7,ij}=-i{\delta }_{i2}{\delta }_{j3}+i{\delta }_{i3}{\delta }_{j2}\\ {\lambda }_{3,ij}\,=\,{{{\rm{diag}}}}\left(1,-1,0\right),\\ {\lambda }_{8,ij}\,=\,\frac{1}{\sqrt{3}}{{{\rm{diag}}}}\left(1,1,-2\right).\end{array}$$

(72)

It is convenient to employ the representation that the spin-1 matrices are purely imaginary, i.e.,

$${({S}^{i})}_{jk}=-i{\epsilon }_{ijk},$$

(73)

which are just the three purely imaginary matrices of Eq. (72),

$${S}^{x}={\lambda }_{7},\quad {S}^{y}=-{\lambda }_{5},\quad {S}^{z}={\lambda }_{2}.$$

(74)

In this representation, S_z is off-diagonal, instead, λ₃ and λ₈ are diagonal. They are spin-quadrupole operators as

$${\lambda }_{3}=-({S}_{x}^{2}-{S}_{y}^{2}),\quad {\lambda }_{8}=-\frac{1}{\sqrt{3}}({S}_{x}^{2}+{S}_{y}^{2}-2{S}_{z}^{2}).$$

(75)

$\frac{1}{2}{\lambda }_{3},\frac{1}{2}{\lambda }_{8}$ span the Catan subalgebra of SU(3). Their common eigenstates are just spin-1 states in the polar bases $\left\vert a\right\rangle$ with a = x, y and z, defined as

$${S}^{a}\left\vert a\right\rangle =0.$$

(76)

Each state can be represented in terms of their eigenvalues of $(\frac{{\lambda }_{3}}{2},\frac{{\lambda }_{8}}{2})$ as

$$\left\vert x\right\rangle :\left(\frac{1}{2},\frac{1}{2\sqrt{3}}\right),\quad \left\vert y\right\rangle :\left(-\frac{1}{2},\frac{1}{2\sqrt{3}}\right),\left\vert z\right\rangle :\left(0,-\frac{1}{\sqrt{3}}\right).$$

(77)

They are vertices of an equilateral triangle in the plane of $(\frac{{\lambda }_{3}}{2},\frac{{\lambda }_{8}}{2})$, which is just the weight diagram of the fundamental representation of the SU(3) group. In other words, $\left\vert x\right\rangle$, $\left\vert y\right\rangle$, and $\left\vert z\right\rangle$ play the role of three colors of quarks. They can represented as [1], i.e., the single box in terms of the Young tableau.

In this representation, Eq. (40) becomes

$$\begin{array}{l}\frac{H}{J}\,=\,\frac{1}{2}\sin \theta \mathop{\sum}\limits_{a=1,3,4,6,8}\mathop{\sum}\limits_{\langle ij\rangle }{\lambda }_{a}(i){\lambda }_{a}(j)\\\qquad\quad+\,\left(\cos \theta -\frac{1}{2}\sin \theta \right)\mathop{\sum}\limits_{a=2,5,7}\mathop{\sum}\limits_{\langle ij\rangle }{\lambda }_{a}(i){\lambda }_{a}(j).\end{array}$$

(78)

At $\theta =\frac{1}{4}\pi$ and $\frac{5}{4}\pi$, this Hamiltonian exhibits an explicit SU(3) symmetry as

$$\frac{H}{J}=\pm \frac{\sqrt{2}}{4}J\mathop{\sum }\limits_{a=1}^{8}\mathop{\sum}\limits_{\langle ij\rangle }{\lambda }_{a}(i){\lambda }_{a}(j),$$

(79)

with each site belonging to the fundamental representation of SU(3). The 6-fold degeneracy comes from the decomposition of,

$$[1]\otimes [1]={[1]}^{* }\oplus [2],$$

(80)

in terms of the Young tableau, where [1]^* is the complex conjugation representation of [1], and [2] is 6D representation of one row with two boxes.

On the other hand, a staggered type of the SU(3) symmetry appears at θ = π/2 and 3π/2. In this case, Eq. (78) becomes

$$\frac{H}{J}=\pm \frac{1}{2}\mathop{\sum}\limits_{\langle ij\rangle }\left(\mathop{\sum}\limits_{a=1,3,4,6,8}{\lambda }_{a}(i){\lambda }_{a}(j)-\mathop{\sum}\limits_{a=2,5,7}{\lambda }_{a}(i){\lambda }_{a}(j)\right),$$

in which the terms of the real and imaginary Gell-mann matrices exhibit the opposite signs. For a bipartite lattice, it can be cast into

$$H/J=\mp \frac{1}{2}\mathop{\sum}\limits_{i\in A,j=i\pm 1\in B}{\lambda }_{a}(i){\left(-{\lambda }_{a}(j)\right)}^{* },$$

(81)

where A and B represent two sublattices. Since $-{\lambda }_{a}^{* }$ is the generators of [1]^*, Eq. (81) is also SU(3) symmetric in which two sublattices lie in the fundamental and anti-fundamental representations, respectively. The 8-fold degeneracy arises from the decomposition of [1] ⊗ [1]^* = • ⊗ [2, 1], where • means the SU(3) singlet and [2, 1] is the 8D adjoint representation of the SU(3) group.

The above two different types of SU(3) symmetries lead to different types of physics. In 1D, the spin-1 chain Eq. (40) with the bilinear and bi-quadratic Heisenberg terms exhibit a variety of quantum phases. Along the line of θ = π/4, each site lies in the fundamental representation of SU(3). The system is a gapless state who critical properties are described by the SU(3)₁ Wess-Zumino-Witten model¹⁶². On the other hand, at θ = − π/2, the system is dimerized breaking the translation symmetry, while two sublattices belong to the fundamental and anti-fundamental representations, respectively. In contrast, at θ = 0, i.e., the spin-1 Heisenberg model with only the bi-linear interaction, the system exhibits the Haldane gap without translation symmetry breaking.

Γ-matrices as quadrupoles of spin-3/2

The spin quadrupole matrices are constructed as follows,

$$\begin{array}{lll}{\Gamma }^{1}=\frac{1}{\sqrt{3}}\left({S}_{x}{S}_{y}+{S}_{y}{S}_{x}\right),\\ {\Gamma }^{2}=\frac{1}{\sqrt{3}}\left({S}_{z}{S}_{x}+{S}_{x}{S}_{z}\right),\\ {\Gamma }^{3}=\frac{1}{\sqrt{3}}\left({S}_{z}{S}_{y}+{S}_{y}{S}_{z}\right),\\ {\Gamma }^{4}={S}_{z}^{2}-\frac{5}{4},\\ {\Gamma }^{5}=\frac{1}{\sqrt{3}}\left({S}_{x}^{2}-{S}_{y}^{2}\right).\end{array}$$

(82)

Remarkably, when S_x,y,z are spin-3/2 matrices, the Γ-matrices defined above anticommute with each other, {Γ^a, Γ^b} = 2δ_ab, forming a representation of the Dirac Γ-matrices. Their commutators are defined as

$${\Gamma }^{ab}=-\frac{i}{2}[{\Gamma }^{a},{\Gamma }^{b}]\quad (1\le a,b\le 5).$$

(83)