text:
\section{Introduction}
\label{sec:intro}
The symmetry principle plays a fundamental role with respect to the laws of nature. It provides an infrastructure and coherence for summarizing physical laws that are independent of any specific dynamics.
Noether's theorem says that every continuous symmetry of the action of a physical system is associated with a corresponding conservation law.
The standard paradigm for describing phase transitions and critical phenomena is Landau's theory of symmetry breaking. The states of matter are classified on the basis of symmetries.
A higher-temperature phase is of a high symmetry characterized by a group $G$, while a lower-temperature phase is of a low symmetry characterized by a subgroup $H\subset G$.
A low-energy effective theory can be constructed in terms of order parameters and is described by all terms that are allowed according to the relevant symmetries.
A hierarchy of symmetries is also widely used in particle physics to understand the dynamics of elementary particles.

Meanwhile, spin-1 quantum magnets are of great interest in physics. One famous example is the Haldane phase in one-dimensional (1D) spin-1 chains\cite{hgap},
in which fractional spin-$1/2$ end states are protected by the spin rotational symmetry in a phenomenon called symmetry-protected topological order\cite{spt1,eespt,spt2,spt3}.
Spin-1 systems are also able to host spin nematic orders in dimensions of $D>1$; such orders are characterized by long-range spin quadrupolar correlations, and the possibility for fractional spinon excitations to coexist with spin nematic orders has also been proposed\cite{s1fermi1}.
Such quantum magnets are widely encountered in various materials, especially transition metal compounds, in which a local $S=1$ magnetic moment can be formed in a cation via Hund's coupling; examples include $3d^{8}$ Ni$^{2+}$ and $3d^{6}$ Fe$^{2+}$.
In this work, we shall reveal several hidden $SU(2)$ symmetries in spin-1 quantum magnets in addition to spin rotational symmetry, and we will study spin-1 quantum systems with the help of the symmetry hierarchy.

For a spin-1 quantum magnet, there are three local states, namely, $|S^{z}=\pm1\rangle$ and $|S^{z}=0\rangle$, and eight independent local Hermitian operators: three spin vector operators, $S^{x}$, $S^{y}$ and $S^{z}$, and five spin quadrupolar operators,
\begin{eqnarray}
Q^{x^{2}-y^{2}} & = & (S^{x})^{2}-(S^{y})^{2},\nonumber \\
Q^{3z^{2}-r^{2}} & = & \frac{1}{\sqrt{3}}\left[2(S^{z})^{2}-(S^{x})^{2}-(S^{y})^{2}\right], \nonumber \\
Q^{xy} & = & S^{x}S^{y}+S^{y}S^{x},  \\
Q^{yz} & = & S^{y}S^{z}+S^{z}S^{y}, \nonumber \\
Q^{zx} & = & S^{z}S^{x}+S^{x}S^{z}. \nonumber
\end{eqnarray}
To illustrate the symmetry hierarchy, we consider a generic two-body interacting Hamiltonian as follows:
\begin{equation}\label{eq:Hgeneric}
\mathcal{H}=\sum_{\langle i,j\rangle}\Big(\sum_{\alpha\beta}J_{\alpha\beta}S^{\alpha}_{i}S^{\beta}_{j}+ \sum_{\mu\nu}J_{\mu\nu}Q^{\mu}_{i}Q^{\nu}_{j}+ \sum_{\alpha\mu}I_{\alpha\mu}S^{\alpha}_{i}Q^{\mu}_{j}\Big),
\end{equation}
where $\langle i,j\rangle$ is a pair of nearest neighboring sites; $\alpha$ and $\beta$ denote $x$, $y$, and $z$; and $\mu$ and $\nu$ denote $x^{2}-y^{2}$, $3z^{2}-r^{2}$, $xy$, $yz$, and $zx$.
The $SO(3)$ spin rotational symmetry is achieved when $I_{\alpha\mu}=0$, $J_{\alpha\beta}=\delta_{\alpha\beta}J_{1}$ and $J_{\mu\nu}=\delta_{\mu\nu}J_{2}$. Furthermore, $\mathcal{H}$ will be $SU(3)$ symmetric when $J_1=J_2$, and the $SU(3)$ group is generated by eight operators $\{\bm{S},\bm{Q}\}$\cite{bbq1d4,SU3ULSL,SU3ULSS}.
The $SO(3)$ model is well studied: a phase diagram consisting of a ferromagnetic phase, a dimerized phase, Haldane phases and a critical phase has been constructed in one dimension\cite{bbq1d1,bbq1d2,bbq1d3,gutzsu3}, and the $SO(3)$ model can host spin nematic ground states in dimensions of $D>1$\cite{s1fermi1,bbq2d1,bbq2d2,bbq2d3,bbq2d4,bbq2d5,bbq2d6,bbq2d7,bbq2d8}.

The model defined in Eq.~\eqref{eq:Hgeneric} has typically been studied in accordance with the $SU(3)\supset SO(3) \cdots$ symmetry hierarchy.
Nevertheless, there are other $SU(2)$ subgroups belonging to the $SU(3)$ group,
and this fact implies the existence of a slice of $SU(2)$ symmetries in addition to the $SO(3)$ spin rotation in spin-1 quantum magnets of which, to the best of our knowledge, the research community is not aware. This situation inspires us to search for Hamiltonians that respect these hidden symmetries;
for this purpose, a new symmetry hierarchy, $SU(3)\supset SU(2)\times U(1)\cdots$, will be adopted to reveal novel states with various low-energy excitations. To describe these states, it is natural to treat all operators $\{\bm{S},\bm{Q}\}$ on an equal footing, which is reminiscent of the ``eight-fold way" in quantum chromodynamics (QCD).

The paper is organized as follows. We reveal hidden $SU(2)$ symmetries  and construct corresponding $SU(2)$ symmetric Hamiltonians in Section~\ref{sec:sym.H}. In the spirit of eight-fold way, we apply flavor-wave mean-field theory to study the $SU(2)\times U(1)$ model and demonstrate that the similar structure in spin spectral functions and spin quadrupole spectral functions serves a clue to the eight-fold way in Section~\ref{sec:MFT}.
In Section~\ref{sec:beyondMFT}, we go beyond the mean-field theory and find an emergent $S=1/2$ gapless quantum spin liquid state coexisting with spin nematic order. We also find Gell-Mann-Okubo formula like gap relations for the hydrodynamic modes in quantum spin-orbial liquid states.
Section~\ref{sec:summary} is devoted to summary.

\section{Hidden $SU(2)$ Symmetries and Hamiltonians}\label{sec:sym.H}

\subsection{Hidden $SU(2)$ Symmetries}
\label{sec:hsu2}
It turns out that there are three hidden $SU(2)$ symmetries in addition to well known spin rotational $SO(3)$ symmetry in a spin-1 system, which are generated as follows (see Appendix.\ref{app:SU2}):
\begin{eqnarray}
SU(2)_{\alpha}&:&\{Q^{zx},S^{y},\frac{1}{2}Q^{x^2-y^2}+\frac{\sqrt{3}}{2}Q^{3z^2-r^2}\},\nonumber\\
SU(2)_{\beta}&:&\{Q^{yz},S^{x},\frac{1}{2}Q^{x^2-y^2}-\frac{\sqrt{3}}{2}Q^{3z^2-r^2}\},\nonumber\\
SU(2)_{\gamma}&:&\{Q^{xy},S^{z},Q^{x^2-y^2}\}.\label{eq:SU2SQ}
\end{eqnarray}
Each set of these generators consists of one component of the spin vector $\bm{S}$ and two components of the spin quadrupole $\bm{Q}$.
Note that these three sets of generators are related to each other by the following cycle: $S^{x}\to S^{y}\to S^{z}\to S^{x}$.
In the remaining part of this work, we shall focus on $SU(2)_{\gamma}$; $SU(2)_{\alpha}$ and $SU(2)_{\beta}$ can then be obtained in accordance with this cycle.

For the $SU(2)_{\gamma}$ symmetry, $S^{z}$ generates spin rotations along the $z$-axis, and the other two generators, $Q^{xy}$ and $Q^{x^2-y^2}$, correspond to {\em two-magnon} processes,
as can be seen from
\begin{subequations}
\begin{eqnarray}
Q^{x^2-y^2}&=&\frac{1}{2}\left[(S^{+})^2+(S^{-})^2 \right],\\
Q^{xy}&=&\frac{1}{2i}\left[(S^{+})^2-(S^{-})^2 \right],
\end{eqnarray}
\end{subequations}
where $S^{\pm}=S^{x}\pm iS^{y}$. Let us define,
\begin{subequations}
\begin{eqnarray}
J^z &=&\frac{1}{2}S^z, \\
J^+ &=&\frac{1}{2}(Q^{x^2-y^2}+iQ^{xy})=\frac{1}{2}(S^+)^2,\\
J^- &=&\frac{1}{2}(Q^{x^2-y^2}-iQ^{xy})=\frac{1}{2}(S^-)^2.
\end{eqnarray}
\end{subequations}
It is easy to verify that $\{J^z,J^{\pm}\}$ satisfy the $SU(2)$ Lie algebra.
Therefore, the spontaneous breaking of the $SU(2)_{\gamma}$ symmetry along the $S^{z}$ direction will give rise to {\em two-magnon} low-energy excitations,
while spontaneous symmetry breaking along the $Q^{xy}$ and $Q^{x^2-y^2}$ directions will give rise to an admixture of {\em one- and two-magnon} excitations, which will tend to restore the $SU(2)_{\gamma}$ symmetry.

The underlying $SU(3)$ structure and the hidden $SU(2)$ symmetries will be more transparent in the Cartesian representation of the spin states: $|x\rangle=i (|1\rangle-|-1\rangle)/\sqrt{2}$, $|y\rangle=(|1\rangle+|-1\rangle)/\sqrt{2}$, and $|z\rangle=-i|0\rangle$.
Then, a spin state can be written as $|\bm{d}\rangle=d^{x}|x\rangle+d^{y}|y\rangle+d^{z}|z\rangle$, where $\bm{d}=\left(d^{x},d^{y},d^{z}\right)$ is a complex vector and the normalization condition is given by $|\bm{d}|^{2}=1$.
The expectation values for $\{\bm{S},\bm{Q}\}$ can be expressed in terms of $\bm{d}$ as follows:
\begin{eqnarray}
\langle S^{\alpha}\rangle & = & -i\epsilon_{\alpha\beta\gamma}\bar{d}^{\beta}d^{\gamma}, \nonumber \\
\langle Q^{\alpha\beta} \rangle|_{\alpha\neq\beta} & = & -(\bar{d}^{\alpha}d^{\beta}+\bar{d}^{\beta}d^{\alpha}), \nonumber \\
\langle Q^{x^{2}-y^{2}}\rangle & = &|d^{y}|^{2}-|d^{x}|^{2}, \nonumber \\
\langle Q^{3z^{2}-r^{2}}\rangle & = & \frac{1}{\sqrt{3}}(2|d^{z}|^{2}-|d^{y}|^{2}-|d^{x}|^{2}),
\end{eqnarray}
where $\bar{d}^{\alpha}$ is the complex conjugate of $d^{\alpha}$ and $\epsilon^{\alpha\beta\gamma}$ is a three-rank antisymmetric tensor. Thus, a spin-1 quantum system can be described by the following path integral:
\begin{equation}\label{eq:pathint}
\mathcal{Z}=\int\mathcal{D}[\bm{d},\bm{\bar{d}}\,] \delta(|\bm{d}|^2-1) e^{-\int_0^{\beta}d\tau \left\{ \sum_{i}\bm{\bar{d}}_i\cdot \partial_{\tau}\bm{d}_i-\mathcal{H} \right\}},
\end{equation}
where the Hamiltonian $\mathcal{H}$ is given by Eq.~\eqref{eq:Hgeneric} with $\{\bm{S},\bm{Q}\}$ replaced with their expectation values.
Now, it is clear that all of the special unitary transformations of $\bm{d}$ give rise to the $SU(3)$ group and that the special unitary transformations of any two components of $\bm{d}$ lead to either $SU(2)_{\alpha}$, $SU(2)_{\beta}$ or $SU(2)_{\gamma}$.

\subsection{$SU(2)$-symmetric Hamiltonians}
\label{sec:su2ham}
Now, we are in a position to construct Hamiltonians in accordance with the $SU(2)_{\gamma}$ symmetry.
A generic spin-1 Hamiltonian can be written in terms of $\{\bm{S},\bm{Q}\}$ in a bilinear form as shown in Eq.~\eqref{eq:Hgeneric}.
Using group theory, one is able to obtain all $SU(2)_{\gamma}$-symmetric two-body interactions (see Appendix.\ref{app:SU2state}).
These $SU(2)_{\gamma}$-symmetric Hamiltonians are linear combinations of the following six terms:
\begin{subequations}\label{eq:H1-6}
\begin{eqnarray}
\mathcal{H}_{1} &=& \sum_{\langle i,j\rangle}S^{z}_{i}S^{z}_{j}+Q^{xy}_{i}Q^{xy}_{j}+Q^{x^{2}-y^{2}}_{i}Q^{x ^{2}-y^{2}}_{j},\\
\mathcal{H}_{2} &=& \sum_{\langle i,j\rangle}S^{x}_{i}S^{x}_{j}+Q^{yz}_{i}Q^{yz}_{j}+S^{y}_{i}S^{y}_{j}+Q^{zx}_{i}Q^{zx}_{j},\\
\mathcal{H}_{3} &=& \sum_{\langle i,j\rangle}Q^{3z^{2}-r^{2}}_{i}Q^{3z^{2}-r^{2}}_{j},\\
\mathcal{H}_{4} &=& \sum_{\langle i,j\rangle}D_{ij}\left[S^{x}_{i}S^{y}_{j}+Q^{yz}_{i}Q^{zx}_{j}-(i\leftrightarrow{}j)\right],\\
\mathcal{H}_{5} &=& \sum_{\langle i,j\rangle}D_{ij}\left[S^{y}_{i}Q^{zx}_{j}+Q^{yz}_{i}S^{x}_{j}-(i\leftrightarrow{}j)\right],\\
\mathcal{H}_{6} &=& \sum_{\langle i,j\rangle}D_{ij}\left[Q^{zx}_{i}S^{x}_{j}+Q^{yz}_{i}S^{y}_{j}-(i\leftrightarrow{}j)\right],
\end{eqnarray}
\end{subequations}
where the $D_{ij}=-D_{ji}=\pm 1$ define a direction along each bond $\langle i,j\rangle$ and are translationally invariant.
The Hamiltonians $\mathcal{H}_{1-6}$ can be classified with respect to the time-reversal symmetry $\mathcal{T}$, the spatial inversion symmetry $\mathcal{I}$, and additional symmetries as summarized in Table \ref{tab:H1-6}.

Further discussions are presented as follows:
(1) $\mathcal{H}_{1}$ consists of the $SU(2)_{\gamma}$ generators $\{S^{z}_i,Q^{xy}_i,Q^{x^{2}-y^{2}}_i\}$ and can be viewed as a spin-1/2 Heisenberg model in the subspace spanned by the local basis $\{|x\rangle_i,|y\rangle_i\}$.
(2) Let us define $\Lambda_8=\sum_{i}Q^{3z^{2}-r^{2}}_{i}$; then, we have $[\Lambda_8,\mathcal{H}_{1}]=[\Lambda_8,\mathcal{H}_{2}]=[\Lambda_8,\mathcal{H}_{3}]=[\Lambda_8,\mathcal{H}_{5}]=0$.
Thus, $\mathcal{H}_{1-3}$ and $\mathcal{H}_{5}$ have an $SU(2)_{\gamma}\times U(1)$ symmetry, where the additional global $U(1)$ symmetry is generated by $\Lambda_8$.
(3) $\mathcal{H}_{1}$ has an additional local $U(1)$ symmetry generated by $Q_{i}^{3z^2-r^2}$.
(4) $\mathcal{H}_{3}$ has an additional local $U(2)$ symmetry that is generated by $\{S^{z}_i,Q^{xy}_i,Q^{x^{2}-y^{2}}_i\}$ and $Q^{3z^{2}-r^{2}}_{i}$.

\begin{table}[htbp]
\caption{Typical $SU(2)_{\gamma}$-symmetric Hamiltonians in Eq.~\eqref{eq:H1-6}, which are classified with respect to time reversal ($\mathcal{T}$), spatial inversion ($\mathcal{I}$), and additional global/local symmetries.}
\label{tab:H1-6}
\renewcommand\arraystretch{1.2}
\begin{tabular}{c|c|c|c|c}
\hline
Hamiltonian &  \hspace{0.5cm}$\mathcal{T}$\hspace{0.5cm} & \hspace{0.5cm}$\mathcal{I}$\hspace{0.5cm} & \hspace{0.6cm}Global\hspace{0.6cm} & \hspace{0.4cm}Local\hspace{0.4cm} \\
\hline
$\mathcal{H}_{1}$ & Yes & Yes & $SU(2)\times U(1)$ & $U(1)$\\
$\mathcal{H}_{2}$ & Yes & Yes & $SU(2)\times U(1)$ & \\
$\mathcal{H}_{3}$ & Yes & Yes & $SU(2)\times U(1)$ & $U(2)$\\
$\mathcal{H}_{4}$ & Yes & No & $SU(2)$ & \\
$\mathcal{H}_{5}$ & No & No & $SU(2)\times U(1)$ & \\
$\mathcal{H}_{6}$ & No & No & $SU(2)$ & \\
\hline
\end{tabular}
\renewcommand\arraystretch{1}
\end{table}

\section{Eight-fold way and flavor-wave mean-field theory: $SU(2)\times U(1)\times\cal{T}\times\cal{I}$ model}\label{sec:MFT}
Regarding the $SU(3)\supset SU(2)\times U(1)$ symmetry hierarchy, we would like to treat all of the $SU(3)$ generators $\{\bm{S},\bm{Q}\}$ on an equal footing.
In this spirit, we study ground states and low-energy excitations with the help of flavor-wave mean-field theory\cite{bbqmodel2,flavor1, flavor2,flavor3,flavor4}, and then discuss how the spectral functions can be exploited to detect the eight-fold way.

First, we minimize the energy functional $\mathcal{H}(\bm{d},\bm{\bar{d}})$ with respect to the complex fields $\bm{d}_i$ to determine the ground state.
Second, we assign three flavors of Schwinger bosons, $a_{n\alpha}(j)$, to each site $j$ on the $n^{th}$ sublattice,
where $\alpha=x,y,z$ refers to the local spin states. For example, $n=1$ for a uniform state, while $n=1,2$ for a bipartite-lattice ordered state.
The operators $\{\bm{S},\bm{Q}\}$ can be written bilinearly in terms of the Schwinger bosons, and the physical Hilbert space can be restored by imposing a single-occupancy condition (see Appendix.\ref{app:FW}).
Third, without loss of generality, we let the Schwinger bosons condense at $a_{n\tilde{x}}$ to obtain ordered states, where $a_{n\tilde{x}}$ and the other two orthogonal components, $a_{n\tilde{y}}$ and $a_{n\tilde{z}}$, are related to $(a_{nx}, a_{ny}, a_{nz})$ by an $SU(3)$ rotation $\Omega_n$ as follows: $(a_{n\tilde{x}}, a_{n\tilde{y}}, a_{n\tilde{z}})^T = \Omega_{n}(a_{nx}, a_{ny}, a_{nz})^T$.
Such an $\Omega_{n}$ is determined by the mean-field vector $\bm{d}$ and enables us to attribute the condensate to $a_{n\tilde{x}}$ alone, while treating $a_{n\tilde{y}}$ and $a_{n\tilde{z}}$ as small fractions. Then, the low-energy Hamiltonian can be bilinearized by the Holstein-Primakoff transformation: $a_{n\tilde{x}}^{\dagger}(j) = a_{n\tilde{x}}(j) = \sqrt{M-a_{n\tilde{y}}^{\dagger}(j)a_{n\tilde{y}}(j)-a_{n\tilde{z}}^{\dagger}(j)a_{n\tilde{z}}(j)}$,
where we will ultimately take $M=1$ for the single-occupancy case.
Expansion in $1/M$ and Bogoliubov transformation will give rise to a diagonalized Hamiltonian in $k$-space (see Appendix.\ref{app:FW}):
$\mathcal{H}=\sum_{m,\bm{k}}\omega_{m}(\bm{k})b^{\dagger}_{m}(\bm{k})b_{m}(\bm{k})+\mathcal{C}$,
where $\omega_{m}(\bm{k})$ is the energy dispersion of the $m$-th flavor-wave branch, $b_{m}(\bm{k})$ is a bosonic Bogoliubov quasiparticle, and $\mathcal{C}$ is a constant.
For a uniform state, $m=1,2$, while for a bipartite-lattice ordered state, $m=1,2,3,4$.
As long as the vector $\bm{d}$ is given by the mean-field theory, we will be able to obtain $\omega_{m}(\bm{k})$ and $b_{m}(\bm{k})$ simultaneously.

In particular, we are interested in Hamiltonians with the time-reversal symmetry $\mathcal{T}$ and the spatial inversion symmetry $\mathcal{I}$, which can be parameterized in terms of three real numbers $K_1$, $K_2$, and $K_3$ as follows:
\begin{equation}\label{eq:H123}
\mathcal{H}=K_{1}\mathcal{H}_{1}+K_{2}\mathcal{H}_{2}+K_{3}\mathcal{H}_{3}.
\end{equation}
Note that the model given in Eq.~\eqref{eq:H123} respects the $SU(2)_{\gamma}\times U(1)$ symmetry rather than the $SU(2)_{\gamma}$ symmetry.
For simplicity, we shall consider bipartite lattices only, including a 1D chain, a square lattice and a cubic lattice.

\subsection{Ground states}
\label{sec:gs}
\begin{figure}[htpb]
\centering
\includegraphics[width=8.4cm]{phasediagram.pdf}
\caption{Mean-field phase diagrams for $SU(2)\times U(1)\times\mathcal{T}\times\mathcal{I}$-symmetric models as defined in Eq.~\eqref{eq:H123} on bipartite lattices.
Two $SU(3)$ points are located at $K_1=K_2=K_3$. (a) Top view of the parameter space ($K_3>0$). (b) Bottom view of the parameter space ($K_3<0$).
Solid lines indicate first-order phase transitions, while dashed lines indicate continuous transitions.
}\label{Fig:phases}
\end{figure}

To explore the ground-state phase diagram, we set $K^{2}_1+K^{2}_2+K^{2}_3=1$, such that the parameter space is a sphere.
Top and bottom views (along the $K_3$ axis) of this sphere are displayed in Fig.~\ref{Fig:phases},
where the mean-field phase diagram is presented.
There are six ordered phases, FQ1, FQ2, FQ3, AFQ1, AFQ2, and AFQ3. Here, FQ refers to a ferro-quadrupolar state, and AFQ refers to an antiferro-quadrupolar state (or, to be precise, a state with a staggered quadrupolar order).
When $K_{1(2,3)}$ is negative and predominates, the ground states are FQ states, while when $K_{1(2,3)}$ is positive and predominates, the ground states are AFQ states.
The solid lines in the phase diagram represent first-order transitions, while the dashed lines represent continuous transitions.
The $SU(3)$ symmetry will be achieved at two points where $K_1=K_2=K_3$. Both $SU(3)$ points are tricritical points. The one with $K_{1,2,3}<0$ corresponds to three phases, FQ1, FQ2 and FQ3,
while the other one, with $K_{1,2,3}>0$, corresponds to AFQ1, AFQ2 and AFQ3.

We use local spin density to illustrate the wavefunctions for various FQ and AFQ states.
If $\bm{d}$ vector is real, the state is time reversal invariant, and the local spin density $|\langle\hat{S}(\hat{n})|\bm{d}\rangle|^{2}$ is invariant under an $SO(2)$ rotation along the axis of $(|d_{x}|, |d_{y}|, |d_{z}|)$, which is so-called $SO(2)$ symmetric pattern.
Otherwise, $|\bm{d}\rangle$ breaks the time-reversal symmetry, $|\langle\hat{S}(\hat{n})|\bm{d}\rangle|^{2}$ will be distorted from an $SO(2)$ symmetric pattern to an $SO(2)$ non-symmetric pattern. Two examples for time reversal invariant state and time reversal breaking state are given in FIG.~\ref{fig:wf0}

\begin{figure}[htpb]
	\includegraphics[width=8.4cm]{distort.pdf}
	\caption{Local spin density for (a) time reversal invariant  and (b) time reversal breaking states. Red bars indicate the directions of $(|d_{x}|, |d_{y}|, |d_{z}|)$. Blue surfaces represent the spin density of a local wavefunction $\Psi$, which is defined as $|\langle\hat{S}(\hat{n})|\Psi\rangle|^{2}$. Here $|\hat{S}(\hat{n})\rangle$ is the spin coherent state pointing along direction $\hat{n}$ and is defined by $(\hat{n}\cdot\bm{S})|\hat{S}(\hat{n})\rangle=|\hat{S}(\hat{n})\rangle$. The local states $|\Psi\rangle$ are chosen for (a) $|\Psi\rangle=(|x\rangle+|y\rangle)/\sqrt{2}$; (b) $|\Psi\rangle=(e^{i\pi/10}|x\rangle+|y\rangle)/\sqrt{2}$. }\label{fig:wf0}
\end{figure}

All the six quadrupolar ordered states are illustrated on square lattice in FIG.~\ref{fig:wf}, where we choose time-reversal invariant ground states to eliminate spin dipolar orders and manifest quadrupolar orders. Notably, dipolar and quadrupolar orders may coexist in a ground state in the FQ1, FQ2, AFQ1 and AFQ2 phases, while only a quadrupolar order exists in the FQ3 and AFQ3 phases.

\begin{figure}[htpb]
	\includegraphics[width=8.4cm]{spinprob.pdf}
	\caption{ Six types of spin quadrupolar orders on square lattice. All the states are time-reversal invariant with real $\bm{d}$ vectors, and red bars indicate the directions of these real $\bm{d}$ vectors. Blue surfaces represent the spin density of a local wavefunction $\Psi$, which is defined as $|\langle\hat{S}(\hat{n})|\Psi\rangle|^{2}$. Here $|\hat{S}(\hat{n})\rangle$ is the spin coherent state pointing along direction $\hat{n}$ and is defined by $(\hat{n}\cdot\bm{S})|\hat{S}(\hat{n})\rangle=|\hat{S}(\hat{n})\rangle$. The local states $|\Psi\rangle$ are chosen for FQ and AFQ phases as follows, (a) FQ1: $|\Psi\rangle=(|x\rangle+|y\rangle)/\sqrt{2}$; (b) AFQ1: $|\Psi_1\rangle=(|x\rangle+|y\rangle)/\sqrt{2}$ and $|\Psi_2\rangle=(|x\rangle-|y\rangle)/\sqrt{2}$; (c) FQ2: $|\Psi\rangle=(|x\rangle+|z\rangle)/\sqrt{2}$; (d) AFQ2: $|\Psi_1\rangle=(|x\rangle+|z\rangle)/\sqrt{2}$ and $|\Psi_2\rangle=(|x\rangle-|z\rangle)/\sqrt{2}$; (e) FQ3: $|\Psi\rangle=|z\rangle$; (f) AFQ3: $|\Psi_1\rangle=(|x\rangle+|y\rangle)/\sqrt{2}$ and $|\Psi_2\rangle=|z\rangle$. Here the subscripts in $\Psi_{1,2}$ refer to sublattices 1 and 2.
	}\label{fig:wf}
\end{figure}

\subsection{Low energy excitations}
\label{sec:exci}
The low-energy excitations can be understood in the framework of the symmetry hierarchy as follows.
(1) The spontaneous symmetry breaking is distinct in the different phases: (a) $SU(2)$ is broken in FQ1 (AFQ1), but $U(1)$ is not (i.e., $SU(2)\times U(1) \to U(1)$); (b) both $SU(2)$ and $U(1)$ are broken in FQ2 (AFQ2) (i.e., $SU(2)\times U(1) \to 1$); and (c) neither $SU(2)$ nor $U(1)$ is broken in FQ3 (AFQ3).
(2) For FQ1, the $\omega^{\mbox{\tiny{FQ1}}}_{1}$ mode is gapless, while the other mode, $\omega^{\mbox{\tiny{FQ1}}}_{2}$, is gapful.
Since $SU(2)$ is broken, the gapless Goldstone mode $\omega^{\mbox{\tiny{FQ1}}}_{1}$ tends to recover the symmetry.
However, $U(1)$ is unbroken, so $\omega^{\mbox{\tiny{FQ1}}}_{2}$ is not required to be gapless as well. The gapless $\omega^{\mbox{\tiny{FQ1}}}_{1}$ mode corresponds to two-magnon excitations, while the gapful $\omega^{\mbox{\tiny{FQ1}}}_{2}$ mode corresponds to one-magnon excitations (see Appendix.\ref{app:FW}).
(3) For FQ2, there are two gapless Goldstone modes, $\omega^{\mbox{\tiny{FQ2}}}_{1}$ and $\omega^{\mbox{\tiny{FQ2}}}_{2}$, because both $SU(2)$ and $U(1)$ are broken. The $\omega^{\mbox{\tiny{FQ2}}}_{1}$ mode is an admixture of one- and two-magnon excitations, while the $\omega^{\mbox{\tiny{FQ2}}}_{2}$ mode consists of one-magnon excitations only.
(4) For FQ3, there are two gapful modes, $\omega^{\mbox{\tiny{FQ3}}}_{1}=\omega^{\mbox{\tiny{FQ3}}}_{2}$, which are related to each other through the $SU(2)$ symmetry. Both of them correspond to one-magnon excitations.
(5) The AFQ1, AFQ2 and AFQ3 phases can be analyzed similarly.

The mean-field ground states and low-energy flavor-wave excitations for these six phases are summarized in Table \ref{tab:summarysix}.

\begin{table*}[htbp]
	\renewcommand\arraystretch{1.3}
	\addtolength{\tabcolsep}{1.0ex}
	\begin{center}
		\caption{Summary of the $SU(2)_{\gamma}\times U(1) \times \mathcal{T} \times \mathcal{I}$-symmetric model defined in Eq.~\eqref{eq:H123}. The parameters $\vartheta$ and $\tilde{\vartheta}$ are given by
$\sin^{2}\vartheta=\frac{2|K_{2}|+K_{1}+K_{3}}{4|K_{2}|+K_{1}+3K_{3}}$ and $\sin^{2}\tilde{\vartheta}=\frac{2K_{2}+K_{1}+K_{3}}{4K_{2}+K_{1}+3K_{3}}$, respectively. $\mathcal{R}(\chi,\theta,\phi,\varphi)=\mbox{diag} (
e^{i\frac{\phi}{2}\sigma^{z}}e^{i\frac{\theta}{2}\sigma^{y}}e^{i\frac{\chi}{2}\sigma^{z}}, e^{i\varphi})$ is an $SU(2)\times U(1)$ rotation. $A_{K}$, $B_{K}$, $C_{K}$, $D_{K}$ and $\gamma(\bm{k})$ are defined as follows:
$A_{K} = \frac{2|K_{2}|-K_{1}+K_{3}-2K_{1}K_{3}/|K_{2}|}{4|K_{2}|+K_{1}+3K_{3}}$,
$B_{K} = \frac{K_{1}+3K_{3}+2K_{3}(K_{1}+K_{3})/|K_{2}|}{4|K_{2}|+K_{1}+3K_{3}}$,
$C_{K} = (3+K_{1}/K_{3})/2$,
$D_{K} = 2K_{2}/K_{3}$, and $\gamma(\bm{k}) = Z^{-1}\sum_{\delta}e^{i\bm{k}\cdot\bm{\delta}}$, where $Z=2D$ is the coordination number and $\bm{\delta}$ is a nearest-neighbor displacement.
``1" refers to one-magnon excitations, ``2" refers to two-magnon excitations, and ``1+2" refers to an admixture of one- and two-magnon excitations.
		 }\label{tab:summarysix}
	\end{center}
	\begin{tabular}{c|c|c|c|c}
		\hline\hline
	       & $\bm{d}$ vector(s)  &  Flavor-wave dispersion & Gap & Magnon \\
		\hline
    FQ1   & $\begin{array}{lll}\bm{d}&=&\mathcal{R}(\chi,\theta,\phi,0)\left(1,0,0\right)^T\end{array}$
    & $\begin{array}{lll}\omega^{\mbox{\tiny{FQ1}}}_{1}(\bm{k}) &=& 2Z|K_{1}|[1-\gamma(\bm{k})]\\
      \omega^{\mbox{\tiny{FQ1}}}_{2}(\bm{k}) &=& Z(|K_{1}|-K_{3})+2ZK_{2}\gamma(\bm{k})\end{array}$
    & $\begin{array}{c}\mbox{gapless}\\ \mbox{gapful}\end{array}$ & $\begin{array}{c}2\\1\end{array}$\\
		\hline
	 FQ2  & $\begin{array}{lll}\bm{d}&=&\mathcal{R}(\chi,\theta,\phi,\varphi)\left(\cos\vartheta,0,\sin\vartheta\right)^T\end{array}$
	 & $\begin{array}{lll} \omega^{\mbox{\tiny{FQ2}}}_{1}(\bm{k}) &=& 2Z|K_{2}|A_{K}(1-\gamma(\bm{k}))\\
      \omega^{\mbox{\tiny{FQ2}}}_{2}(\bm{k}) &=& 2Z|K_{2}|\sqrt{[1-\gamma(\bm{k})][1+B_{K}\gamma(\bm{k})]} \end{array}$
    & gapless & $\begin{array}{c}1+2\\1\end{array}$\\
	  \hline
	 FQ3  & $\begin{array}{lll}\bm{d}&=&\mathcal{R}(0,0,0,\varphi_{i})\left(0,0,1\right)^T\end{array}$
	 & $\begin{array}{lllll}\omega^{\mbox{\tiny{FQ3}}}_{1}(\bm{k}) & = & \omega^{\mbox{\tiny{FQ3}}}_{2}(\bm{k}) & = & 2Z[|K_{3}|+K_{2}\gamma(\bm{k})] \end{array}$ & gapful & $1$\\
		\hline
    AFQ1 & $\begin{array}{lll}\bm{d}_{1} & = & \mathcal{R}(\chi,\theta,\phi,0)\left(1,0,0\right)^T \\ \bm{d}_{2} & = & \mathcal{R}(\chi,\theta,\phi,0)\left(0,1,0\right)^T \end{array}$
	  & $\begin{array}{lllll} \omega^{\mbox{\tiny{AFQ1}}}_{1}(\bm{k}) & = & \omega^{\mbox{\tiny{AFQ1}}}_{2}(\bm{k}) & = & 2ZK_{1}\sqrt{1-\gamma^{2}(\bm{k})}\\
      \omega^{\mbox{\tiny{AFQ1}}}_{3}(\bm{k}) & = & \omega^{\mbox{\tiny{AFQ1}}}_{4}(\bm{k}) & = & Z(K_{1}-K_{3})\end{array}$
     & $\begin{array}{c}\mbox{gapless}\\ \mbox{gapful}\end{array}$ & $\begin{array}{c}2\\1\end{array}$\\
		\hline
		AFQ2  & $\begin{array}{lll}\bm{d}_{1} & = & \mathcal{R}(\chi,\theta,\phi,\varphi)\left(\cos\tilde{\vartheta},0,\sin\tilde{\vartheta}\right)^T \\ \bm{d}_{2} & = & \mathcal{R}(\chi,\theta,\phi,\varphi)\left(\cos\tilde{\vartheta},0,-\sin\tilde{\vartheta}\right)^T \end{array}$
		& $\begin{array}{lll} \omega_{1}^{\mbox{\tiny{AFQ2}}} &=& 2ZK_{2}A_{K}[1-\gamma(\bm{k})]\\
      \omega_{2}^{\mbox{\tiny{AFQ2}}} &=& 2ZK_{2}A_{K}[1+\gamma(\bm{k})]\\
      \omega_{3}^{\mbox{\tiny{AFQ2}}} &=& 2ZK_{2}\sqrt{[1+\gamma(\bm{k})][1-B_{K}\gamma(\bm{k})]}\\
      \omega_{4}^{\mbox{\tiny{AFQ2}}} &=& 2ZK_{2}\sqrt{[1-\gamma(\bm{k})][1+B_{K}\gamma(\bm{k})]}\end{array}$
     & gapless & $\begin{array}{c}1+2\\1+2\\1\\1\end{array}$\\
		\hline
		AFQ3  & $\begin{array}{lll}\bm{d}_{1} & = & \mathcal{R}(0,0,0,\varphi_{i})\left(0,0,1\right)^T \\ \bm{d}_{2} & = & \mathcal{R}(\chi_{i},\theta_{i},\phi_{i},0)\left(1,0,0\right)^T \end{array}$
		& $\begin{array}{lll} \omega^{\mbox{\tiny{AFQ3}}}_{1}(\bm{k}) &=& Z(K_{3}-K_{1})\\
      \omega^{\mbox{\tiny{AFQ3}}}_{2}(\bm{k}) &=& 0+O(\bm{k}^{3})\\
      \omega^{\mbox{\tiny{AFQ3}}}_{3}(\bm{k}) &=& Z[K_{3}\sqrt{C_{K}^{2}-D_{K}^{2}\gamma^{2}(\bm{k})}+\frac{K_{1}-K_{3}}{2}]\\
      \omega^{\mbox{\tiny{AFQ3}}}_{4}(\bm{k}) &=& Z[K_{3}\sqrt{C_{K}^{2}-D_{K}^{2}\gamma^{2}(\bm{k})}+\frac{K_{3}-K_{1}}{2}]\end{array}$ & $\begin{array}{c}\mbox{gapful}\\ \mbox{gapless} \\ \mbox{gapful} \\ \mbox{gapful}\end{array}$ & $\begin{array}{c}2\\2\\1\\1\end{array}$ \\
		\hline\hline
	\end{tabular}
\end{table*}

\subsection{Spectral functions: A clue to the eight-fold way}
\label{sec:spec}
Inelastic neutron scattering measures the spin spectral function in $(\bm{q},\omega)$ space, which is defined as $$S^{\alpha\beta}(\bm{q},\omega)=\mbox{Im}\left\{i\int_{0}^{\infty}dt{}e^{i\omega{}t}\langle[S ^{\alpha}(\bm{q},t),S^{\beta}(-\bm{q},0)]\rangle\right\},$$
where we have set $\mu{}g_{B}=1$ for simplicity.
At zero temperature, $S^{\alpha\beta}(\bm{q},\omega)$ depends on the choice of the ground state. However, the spectral function, $$S(\bm{q},\omega)\equiv S^{xx}(\bm{q},\omega)+S^{yy}(\bm{q},\omega)+S^{zz}(\bm{q},\omega),$$
does not change qualitatively within a single phase.
On the other hand, resonant inelastic X-ray scattering (RIXS) measures two-magnon processes, which is described by spin quadrupole spectral functions, $$Q^{\mu\nu}(\bm{q},\omega)=\mbox{Im}\left\{i\int_{0}^{\infty}dt{}e^{i\omega{}t}\langle[Q ^{\mu}(\bm{q},t),Q^{\nu}(-\bm{q},0)]\rangle\right\},$$ and $\mu$ and $\nu$ denote $xy,\ yz,\ zx,\ x^{2}-y^{2}$ and $3z^{2}-r^{2}$. Similar to $S(\bm{q},\omega)$, the spin quadrupole spectral function, $$Q(\bm{q},\omega)\equiv \sum_{\mu}Q^{\mu\mu}(\bm{q}, \omega),$$ does not change qualitatively within a single phase as well.
Therefore, the spin spectral function $S(\bm{q},\omega)$ and the spin quadrupole spectral function $Q(\bm{q},\omega)$ can be used to detect flavor waves and distinguish the various FQ and AFQ phases.

In the flavor-wave mean-field theory, all the degenerate ground states can be obtained from one of them by an $SU(2)\times U(1)$ rotation of the $\bm{d}$-vector. We parameterize a general $SU(2)$ rotational matrix as
\begin{equation}\label{eq:R_ri}
R = r_{0}\sigma^{0} + i\sum_{n=1}^{3}r_{n}\sigma^{n},
\end{equation}
where $\sigma^{0}$ are identity as well as $\sigma^{1}$, $\sigma^{2}$ and $\sigma^{3}$ are three Pauli matrices. Here $\hat{r}=\{r_{0},r_{1},r_{2},r_{3}\}$ is a 4-dimensional real vector with $\sum_{n=0}^{3}r^{2}=1$. Thus, apart from a global phasefactor, two $\bm{d}$ vectors of two degenerate ground states, $\bm{d}$ and $\bm{d}^{\prime}$, are related by a $3\times 3$ matrix as follows,
\begin{equation}
\bm{d}^{\prime} = \left(\begin{array}{cc}
e^{-i\phi}R & 0_{2,1}\\
0_{1,2} & 1 \\
\end{array}\right)\bm{d},
\end{equation}
where $0_{1,2}$($0_{2,1}$) is a $1\times{}2$($2\times{}1$) zero matrix. In the $SU(3)$ Schwinger boson representation, the expression for the spin operators $\bm{S}$ depends on the parameters $r_{0-3}$

Thus these spectral functions can be evaluated for each FQ or AFQ state; these functions are distinct in different phases but do not qualitatively change within a single phase.
Moreover, $S(\bm{q},\omega)$ and $Q(\bm{q},\omega)$ {\em share the same structure} as long as the elementary excitations are flavor waves, as demonstrated in Fig.~\ref{fig:Sqqw}.
Namely, $Q(\bm{q},\omega)$ has the same dispersion as $S(\bm{q},\omega)$, and difference between them is in the spectral weight. {\em This similarity provides evidence for the underlying $SU(3)$ structure and serves as a clue to the eight-fold way.}
The details of these spectral functions for all FQ and AFQ phases can be found in see Appendix.\ref{app:sf}.

\begin{figure*}[htbp]
	\centering
	\includegraphics[scale=0.62]{SQqw2s.pdf
	\caption{The spin spectral functions $S(\bm{q},\omega)$ and the spin quadrupole spectral functions $Q(\bm{q},\omega)$ for the FQ and AFQ phases. Here, we set $K_{2}=K_{3}=0.1K_{1}$ for FQ1 and AFQ1, $K_{1}=K_{3}=0.1K_{2}$ for FQ2 and AFQ2 and $K_{1}=K_{2}=0.1K_{3}$ for FQ3 and AFQ3.}\label{fig:Sqqw}
\end{figure*}

\section{Beyond the mean-field theory}\label{sec:beyondMFT}
\subsection{Effective Hamiltonian in the AFQ3 phase: A possible emergent $S=1/2$ gapless spin liquid}
\label{sec:effqsl}
In the mean-field solution, the AFQ3 ground states are locally degenerate inside a bulk energy gap.
This huge degeneracy arises from the unperturbative Hamiltonian $K_3\mathcal{H}_{3}$, of which ground state subspace is spanned by the local basis $\{|x\rangle_i,|y\rangle_i\}$ on sublattice-2 and $\{|z\rangle_i\}$ on sublattice-1, as shown in FIG.\ref{fig:heff}(a).
And the degeneracy is expected to be lifted by a small but finite $K_1$ and $K_2$.
To address this case and go beyond the mean-field theory, we consider perturbations of up to the third order in the limit $K_3\gg |K_{1(2)}|$, where the unperturbative energy gap is about $ZK_3$.
What we need is to include all the possible perturbations of $\mathcal{H}_1$ and $\mathcal{H}_2$ and project the states into the subspace of unperturbative ground states by a projector $\mathcal{P}$.

We begin with two-site Hamitonians $h_{i}$, $i=1,2,3$, which is the simplest case of the Hamiltonians $\mathcal{H}_{i}$ defined in Eq.~\eqref{eq:H1-6}, and consider their matrix elements. Firstly, $h_{3} = Q^{3z^{2}-r^{2}}_{1}Q^{3z^{2}-r^{2}}_{2}$, which is diagonal in the basis $\{|\sigma_1\sigma_2\rangle\}$, where $\sigma_{1,2}=x,y,z$.
Secondly, $h_{1}=S^{z}_{1}S^{z}_{2}+Q^{xy}_{1}Q^{xy}_{2}+Q^{x^{2}-y^{2}}_{1}Q^{x^{2}-y^{2}}_{2}$, and we have
\begin{equation}\label{eq:h1}
h_{1}|\sigma_{1}\sigma_{2}\rangle = \left\lbrace \begin{array}{cl}
|\sigma_{1}\sigma_{2}\rangle, & \sigma_{1}=\sigma_{2} \neq{}z, \\
2|\sigma_{2}\sigma_{1}\rangle, & \sigma_{1}=x(y),\sigma_{2}=y(x), \\
0, & \mbox{otherwise}.
\end{array}\right.
\end{equation}
Thirdly, $h_{2}={}S^{x}_{1}S^{x}_{2}+S^{y}_{1}S^{y}_{2}+Q^{yz}_{1}Q^{yz}_{2}+Q^{zx}_{1}Q^{zx}_{2}$, we have
\begin{equation}\label{eq:h2}
h_{2}|\sigma_{1}\sigma_{2}\rangle = \left\lbrace \begin{array}{cl}
2|\sigma_{2}\sigma_{1}\rangle, & \sigma_{1} =x,y(z)\  \mbox{and}\  \sigma_{2}=z(x,y), \\
0, & \mbox{otherwise}.
\end{array}\right.
\end{equation}
Taking into account all the nonzero matrix elements, the leading perturbation is of the third order and the effective Hamiltonian reads
\begin{equation}
\mathcal{H}_{eff} =	\frac{K_{1}K^{2}_{2}}{K_{3}^{2}} \mathcal{P}\mathcal{H}_{2}\mathcal{H}_{1}\mathcal{H}_{2}\mathcal{P}.
\end{equation}
Regarding the projector $\mathcal{P}$, $\mathcal{H}_{eff}$ can be written as a $\mathcal{H}_{1}$ type Hamiltonian defined on one of the sublattices.

As an example, consider a square lattice; the spins have a quadrupolar order on one of the two sublattices, and we have the following effective Hamiltonian on the other sublattice:
\begin{eqnarray}
\mathcal{H}_{eff} & = & \tilde{J}_1 \sum_{\langle ij\rangle_1} \mathcal{P}\left(S^{z}_{i}S^{z}_{j}+Q^{xy}_{i}Q^{xy}_{j}+Q^{x^{2}-y^{2}}_{i}Q^{x^{2}-y^{2}}_{j}\right)\mathcal{P} \nonumber \\
&& + \tilde{J}_2 \sum_{\langle ij\rangle_2}\mathcal{P}\left(S^{z}_{i}S^{z}_{j}+Q^{xy}_{i}Q^{xy}_{j}+Q^{x^{2}-y^{2}}_{i}Q^{x^{2}-y^{2}}_{j}\right)\mathcal{P}\nonumber,
\end{eqnarray}
where $\langle ij\rangle_{1(2)}$ denotes a pair of (next) nearest neighboring sites on the sublattice, $\tilde{J}_1=2K_{1}K^{2}_{2}/K_{3}^{2}$ and $\tilde{J}_2=K_{1}K^{2}_{2}/K_{3}^{2}$.
Note that the relation $\tilde{J}_2/\tilde{J}_1=1/2$ is guaranteed by the fact that there are two paths contribute to $\tilde{J}_1$ while only one path contributes to $\tilde{J}_2$, as illustrated in FIG.~\ref{fig:heff}

\begin{figure}[htpb]
	\centering
	\includegraphics[scale=0.23]{hsu2afq3.pdf}
	\caption{Solid circles form sublattice-$1$ where all the spins are in the $|z\rangle$ state and open circles form sublattice-$2$ where all the spins are in $|x\rangle$ or $|y\rangle$ state. (a) Emergent $\tilde{J}_{1}-\tilde{J}_{2}$ super square lattice with black solid lines for $\tilde{J}_{1}$ bonds and blue curve lines for $\tilde{J}_{2}$ bonds. (b) one path to achieve $J_{2}$ bond. (c) two paths to achieve $\tilde{J}_{1}$ bonds. (b) and (c) indicate that $\tilde{J}_{2}/\tilde{J}_{1}=0.5$.
	}\label{fig:heff}
\end{figure}

Thus, there emerges an effective spin-$1/2$ $\tilde{J}_1$-$\tilde{J}_2$ Heisenberg model constructed by the $SU(2)_{\gamma}$ generators.
When $K_1<0$, $\tilde{J}_{1(2)}<0$, and the ground state is of ferromagnetic order. When $K_1>0$, $\tilde{J}_1=2\tilde{J}_2>0$, a gapless quantum spin liquid (QSL) ground state~\cite{RMPQSL} may be favored, as suggested by variational Monte Carlo~\cite{j1j2_gapless_vmc, j1j2_gapless_vmc2, j1j2_gapless_vmc3} and density matrix renormalization group (DMRG)~\cite{j1j2_gapless_dmrg}. It is worth noting that the ground states of spin-$1/2$ $J_1$-$J_2$ Heisenberg model still remains a controversial issue, and various numerical approaches lead to several contradictory results. Other possible ground states include N\'{e}el order\cite{j1j2_neel_peps, j1j2_neel_peps2} (by PEPS), gapful spin liquid~\cite{j1j2_gapped_dmrg} (by DMRG), and plaquette valence-bond solid.\cite{j1j2_PVB_ed, j1j2_PVB_ed2, j1j2_PVB_trg, j1j2_PVB_dmrg} (by exact diagonalization, tensor network, and DMRG). Despite of the controversy, we would like to propose that a QSL state with effective spin $S=1/2$ (or other ordered states) may coexist with gapful quadrupolar order in a quantum spin-1 system.

\subsection{Hydrodynamics modes in a quantum spin-orbital liquid: In analogy to QCD}
\label{sec:hydro}
Now we consider a quantum spin-orbital liquid ground state in the $SU(2)\times U(1)$ symmetric model, where neither $SU(2)$ nor $U(1)$ symmetry is broken and all the correlation functions of the operators $\{\bm{S},\bm{Q}\}$ are short ranged. If there does exist such a spin-orbital liquid state, the low energy excitations can be described by the fields of $\{\bm{S},\bm{Q}\}$ and can be classified according to the $SU(2)\times U(1)$ symmetry. Because of the short ranged correlation, these excitations are gapful. In analogy to QCD, we will demonstrate below that these excitation gaps must satisfy some relations due to the symmetry hierarchy $SU(3)\supset SU(2)\times U(1)$.

An interesting observation is that the three local spin states $|x\rangle,|y\rangle,|z\rangle$ can be naturally in analogy with the $u,d,s$ quarks in particle physics in the fundamental representation of $SU(3)$ symmetry as follows,
\begin{equation}\label{eq:uds}
\renewcommand\arraystretch{1.25}
\begin{array}{lll}
|x\rangle & \longleftrightarrow & u, \\
|y\rangle & \longleftrightarrow & d, \\
|z\rangle & \longleftrightarrow & s.
\end{array}
\end{equation}
The quark model is a successful theory of the strong interaction, which is known as QCD. According to Gell-Mann's argument: (1) There exists an additive quantum number called strangeness is conserved in addition to isospin $SU(2)$ symmetry; (2) in very strong interactions region, the symmetry is $SU(3)$ rather then the $SU(2)\times{}U(1)$; (3) in medium strong interactions region, the $SU(3)$ breaks into $SU(2)\times{}U(1)$, i.e., isospin $SU(2)$ and hypercharge $U(1)$. This symmetry hierarchy is exactly as the same as what we discuss in our spin-1 systems.
Moreover, the $SU(3)$ octet $\{\bm{S},\bm{Q}\}$ can be mapped to high energy particles, e.g., the light spin-0 mesons, in addition to the $SU(3)$ triplet mapping in Eq.~\eqref{eq:uds} as follows,
\begin{equation}\label{eq:mesons}
\renewcommand\arraystretch{1.25}
\begin{array}{llll}
\frac{1}{2}(iS^{x}-Q^{yz})& \longleftrightarrow & K^{0} &  (d\bar{s}), \\
\frac{1}{2}(-iS^{y}-Q^{zx})& \longleftrightarrow & K^{+} & (u\bar{s}), \\
\frac{1}{2}(iS^{x}+Q^{yz})& \longleftrightarrow & \bar{K}^{0} & (s\bar{d}), \\
\frac{1}{2}(-iS^{y}+Q^{zx})& \longleftrightarrow & K^{-} & (s\bar{u}), \\
\frac{1}{2}(iS^{z}-Q^{xy}) & \longleftrightarrow & \pi^{+} & (u\bar{d}), \\
\frac{1}{2}(-iS^{z}-Q^{xy}) & \longleftrightarrow & \pi^{-} & (d\bar{u}), \\
Q^{x^{2}-y^{2}}  & \longleftrightarrow & \pi^{0} & (u\bar{u},d\bar{d}), \\
Q^{3z^2-r^2}  & \longleftrightarrow & \eta & (u\bar{u},d\bar{d},s\bar{s}).
\end{array}
\end{equation}

For the shorthand, let us define the following operators,
\begin{equation}\label{eq:hmodes}
\renewcommand\arraystretch{1.25}
\begin{array}{lll}
K^{0} &=& (iS^{x}-Q^{yz})/2,\\
K^{+} &=& (-iS^{y}-Q^{zx}),\\
\bar{K}^{0} &=& (-iS^{x}-Q^{yz})/2,\\
K^{-} &=& (iS^{y}-Q^{zx})/2,\\
\pi^{+} &=& (iS^{z}-Q^{xy})/2,\\
\pi^{-} &=& (-iS^{z}-Q^{xy})/2,\\
\pi^{0} &=& Q^{x^{2}-y^{2}},\\
\eta &=& Q^{3z^2-r^2}.
\end{array}
\end{equation}
Thus, in the hydrodynamic limit, $\bm{q}\to 0$ and $\omega \to 0$, the collective modes in a quantum spin-orbital liquid can be described by the fields of $K^{0}$, $\bar{K}^{0}$, $K^{\pm}$, $\pi^{0}$, $\pi^{\pm}$ and $\eta$.
And we expect that each hydrodynamic mode has an excitation gap $M$. The excitaiton gap is nothing but the mass of the corresponding particle.
By the $SU(2)\times U(1)$ symmetry, we deduce that these gaps must satisfy the relations,
\begin{subequations}
\begin{eqnarray}
M_{K^{0}} & = & M_{K^{+}}, \\
M_{\bar{K}^{0}} & = & M_{K^{-}}, \\
M_{\pi^{\pm}} & = & M_{\pi^{0}} = M_{\pi},
\end{eqnarray}
as well as the famous Gell-Mann-Okubo formula \cite{Georgi},
\begin{equation}
2(M_{K^{+}}+M_{K^{-}})=3M_{\eta}+M_{\pi}.
\end{equation}
\end{subequations}

\section{Summary}
\label{sec:summary}
In summary, we have revealed hidden $SU(2)$ symmetries in spin-1 quantum magnets, studied them in accordance with the $SU(3)\supset SU(2)\times U(1)$ symmetry hierarchy, demonstrated novel emergent phenomena, and found some clues to the emergent eight-fold way. These $SU(2)$ symmetries may be realized in cold atoms as well as $d^8$ and/or $d^6$ electrons with the proper specific choices of the spin-orbital couplings.

\section{Acknowledgement}
We would like to thank Dong-Hui Xu, Hong-Hao Tu, Gang Chen, Hong Yao and Zheng-Yu Weng for helpful discussions.
This work is supported in part by National Key Research and Development Program of China (No.2016YFA0300202), National Natural Science Foundation of China (No.11774306), the Strategic Priority Research Program of Chinese Academy of Sciences (No. XDB28000000) and the Fundamental Research Funds for the Central Universities in China.
JJM. is supported by China Postdoctoral Science Foundation(Grant No.2017M620880) and the National NaturalScience Foundation of China (Grant No.1184700424).


text:
\section{Acknowledgements}
We are very appreciative of Haim Sompolinsky for his advice on this work. We thank Sean Escola, Rainer Engelken, and Samuel Muscinelli for comments on the manuscript. We thank Sean Escola and James Murray for early conversations. The authors were supported by the Gatsby Charitable Foundation and NSF NeuroNex award DBI-1707398. A. L.-K. was also supported by the McKnight Endowment Fund.

redpajama_set_name: RedPajamaArXiv


text:
\section{Introduction}
Flat electronic bands are susceptible to a variety of instabilities driven by electron-electron correlations and band filling.  In particular, nearest-neighbor hopping on a two-dimensional (2D) kagome lattice guarantees flat bands while also hosting Dirac band crossings and band touchings that impart a non-trivial topology.  Recent discoveries of superconductivity \cite{Ortiz20} and charge ordering \cite{Jiang21} in $A$V$_3$Sb$_5$, and itinerant ferromagnetism (FM) with large (topological) anomalous Hall response in Co$_3$Sn$_2$S$_2$ \cite{Liu18} have elevated interest in kagome metals as an adaptable system to study the interplay of topology, superconductivity, magnetism, and other charge instabilities \cite{Yin22}.

$R$Mn$_6$Sn$_6$ (where $R$ is a rare-earth) materials comprise an interesting class of magnetic kagome metals \cite{Ghimire20,Yin20,Ma21,Li21,Dhakal21,Kabir22}.  The metallic Mn kagome layers have robust itinerant FM order that can be manipulated by interleaved FM $R$ triangular layers through $R$--Mn antiferromagnetic coupling.  In TbMn$_6$Sn$_6$, $R$--Mn coupling and the uniaxial anisotropy of the Tb ions forces Mn moments to orient perpendicular to the kagome layer \cite{Idrissi91}, creating an ideal scenario for a Chern insulator where spin-orbit coupling (SOC) gaps out spin-polarized Dirac band crossings \cite{Yin20}.  In $R$Mn$_6$Sn$_6$ and other kagome metals with FM layers, such as FeSn \cite{Kakihana18, Kang20, Ye20, Han21} and Fe$_3$Sn$_2$ \cite{Kida11,Lin18,Ye19}, the connection between itinerant FM order and flat electronic bands, evidence of the magnetic coupling to Dirac fermions, and the role of magnetic fluctuations in topological phenomena are all open questions.

These key questions can be addressed using inelastic neutron scattering (INS) measurements to probe magnetic excitations and determine their coupling to kagome electronic bands. The high energy of flat and optical magnon bands in FM kagome lattices, which also host topological magnonic features \cite{Onose10,Mook14,Mook14_2}, make them susceptible to Landau damping by particle-hole excitations (Stoner continuum).  The heavy damping that obscures the observation of these modes in a variety of FM kagome metals supports this hypothesis \cite{Zhang20,Xie21,Do21,Riberolles22}.

Here, we use INS measurements on TbMn$_6$Sn$_6$ to provide key insights into the high energy magnetic excitations of a FM kagome metal.  Whereas the low-energy acoustic magnon modes are well-defined collective excitations, we observe unusual excitations at high energies that are described by short-ranged spin correlations on a hexagonal plaquette in the kagome layer.  Excitations at the K-point of the Brillouin zone consist of spin correlations on a hexagonal plaquette that are derived from the expected Wannier states associated with a flat magnon mode \cite{Bergman08}.  Excitations at the zone center ($\Gamma'$) exhibit chiral antiferromagnetic correlations on the hexagonal plaquette that are commonly associated with magnetic frustration \cite{Inami00, Grohol05}.  Density-functional theory (DFT) calculations assess that the overdamped character of these modes can originate from Landau damping \cite{Fawcett88, Diallo09, Chen20} which is not likely to arise from massive Dirac fermions lying close to the Fermi energy \cite{Yin20,Ma21,Li21,Liu21,Gu22}.  This discovery of chiral magnetic cluster excitations raises new questions about their itinerant character and their coupling to other novel electronic features of the kagome lattice, such as orbital plaquette currents \cite{Mielke22}.

 \section{Experimental Details}
Single crystals of Tb166 were grown from excess Sn using the flux method as previously described \cite{Riberolles22}. INS measurements were performed on the Wide Angular-Range Chopper Spectrometer (ARCS) at the Spallation Neutron Source at Oak Ridge National Laboratory. An array of nine crystals with a total mass of 2.56 grams was co-aligned with the ($H$,0,$L$) scattering plane set horizontally, and attached to the cold head of a closed-cycle-refrigerator. The data were collected at the base temperature of 5 K using incident energies of $E_i =$ 250 and 500 meV.  For each $E_i$ measurement, the sample was rotated around the vertical axis to increase the \textbf{q} coverage. The neutron scattering data are described using the momentum transfer in hexagonal reciprocal lattice units, ${\bf q}(H,K,L) = \frac{2\pi}{a}\frac{2}{\sqrt{3}}(H\hat{a}^*+K\hat{b}^*)+\frac{2\pi}{c}L\hat{z}$. The INS data are presented in terms of the orthogonal vectors $(1,0,0)$ and $(-1,2,0)$, as shown in Fig.~\ref{dispersion}(b). We describe the data with reference to special points in the 2D Brillouin zone; $\Gamma$--(0,0), M--($\frac{1}{2}$,0), and K--($\frac{1}{3}$,$\frac{1}{3}$). The INS data are displayed as intensities that are proportional to the spin-spin correlation function $S({\bf q},E)$, where $E$ is the energy. To improve statistics, the data have been symmetrized with respect to the crystallographic space group P6/$mmm$.

\begin{figure}
\includegraphics[width=1.0\linewidth]{dispersion_multi_panel_v2.png}
\caption{\footnotesize (a) Ferrimagnetic structure of Tb166 with key interlayer magnetic interactions indicated by black arrows.  (b) 2D hexagonal Brillouin zone showing reciprocal lattice vectors $\textbf{a}^{*}=(1,0)$ and $\textbf{b}^{*}=(0,1)$ and representative high symmetry $\Gamma$ (black circles), $\Gamma'$ (empty circles), M (blue circles), and K (gold circles) points. INS data are described with the orthogonal vectors (1,0) and (-1,2). Red arrows correspond to the three reciprocal space slices shown in panels (c)--(h). Panels (c)--(e) show the intensity of slices through the INS data along the $(1-K,2K)$ (slice \#1), $(H,0)$ (slice \#2), and $(H-0.5,1)$ (slice \#3) directions, respectively, after averaging over $L=0-7$. In panels (c) and (d), data with $E<150$ meV ($> 150$ meV) were collected with $E_i=250$ meV (500 meV), respectively. Data in (e) were collected with $E_i=500$ meV. Panels (f)-(h) show identical slices as (c)--(e) calculated from linear spin wave theory using the parameters described in the main text.  In panels (c)--(g), the solid and dashed pink lines correspond to model dispersions with $L=0$ and $L=1/2$, respectively.  In panel (h), blue, pink, green, and white lines label the model dispersions of the Tb, AE/AO, OE/OO, and FE/FO modes, respectively, as described in the text.  White (yellow) dashed rectangles outline overdamped cluster excitations at the $\Gamma'$~(K) point, respectively.}
\label{dispersion}
\end{figure}

\section{INS data}
The magnetic sublattice of TbMn$_6$Sn$_6$ consists of stacked Mn kagome and Tb triangular layers, as shown in Fig.~\ref{dispersion}(a).  The FM Mn and Tb layers couple antiferromagnetically, resulting in a ferrimagnetic ground state with moments pointing along the $c$-axis.  Within linear spin wave theory, TbMn$_6$Sn$_6$ has seven spin wave branches. We label these branches as acoustic--even (AE), acoustic--odd (AO), optical--even (OE), optical--odd (OO), flat--even (FE), flat--odd (FO), and Tb.  As there are two kagome layers in the unit cell, the even (odd) branches correspond to in-phase (out-of-phase) precession of Mn moments in adjacent kagome layers and have strong neutron intensity in Brillouin zones with $L=even$ ($L=odd$), respectively. The acoustic and optical branches are found below and above the K-point Dirac magnon crossing, respectively, and represent in-phase and out-of-phase precession of the three Mn moments within the unit cell of a single kagome layer. AE and AO branches have strong neutron intensity in $\Gamma$ zones with $H=even$ and $K=even$, whereas OE, OO, FE, and FO branches are strongest in $\Gamma'$ zones with $H=odd$ and $K=odd$ or $H+K=odd$, as shown in Fig.~\ref{dispersion}(b).

\begin{figure*}
\includegraphics[width=1.\linewidth]{constE_overlay.png}
\caption{\footnotesize Constant energy slices of the TbMn$_6$Sn$_6$ spin excitations over incremental energy ranges from (a) 50--70 meV, (b) 70--90 meV, (c) 90--110 meV, (d) 110--130 meV, (e) 130--150 meV, (f) 150--170 meV, (g) 170--190 meV, (h) 190--210 meV, (i) 210--230 meV, (j) 230--250 meV.  Data were collected at $E_i=250$ meV for panels (a)--(e) and $E_i=500$ meV for panels (f)--(j). All cuts are averaged over an $L$-range from -7 to 7 rlu. Hexagonal Brillouin zone boundaries are shown and $\Gamma$ and $\Gamma'$ zone centers are indicated by filled and empty circles, respectively.}
\label{constE}
\end{figure*}

Previous INS experiments were conducted on smaller samples and mapped out the lowest-lying AE, AO and Tb branches below 125 meV \cite{Riberolles22}.  These spin waves possess sharp, dispersive excitations throughout the Brillouin zone, indicative of their collective nature, and are well represented by a Heisenberg model consisting of intralayer and interlayer pairwise exchange interactions and single-ion anisotropy terms, as described in the Supplementary Material (SM) \cite{SM}.

In Ref.~\cite{Riberolles22}, we were unable to clearly observe the OE, OO, FE, and FO branches due to the small sample volume and increasingly broad line shapes encountered at higher energies. The measurements reported here were performed with a larger sample volume and higher incident energies, revealing significant magnetic spectral weight up to 250 meV that accounts for the missing branches in the previous data. However, as previous reports hinted, these higher energy features are {\it incoherent} (broad in both momentum and energy), unlike the collective nature of the AO, AE, and Tb branches.

Figures \ref{dispersion}(c)--(e) show slices of the data along different reciprocal space directions within the kagome layers, as indicated in Fig.~\ref{dispersion}(b). The data are averaged over $L$ to improve statistics, resulting in the simultaneous observation of even and odd branches. The data are compared to the Heisenberg model dispersions shown as pink lines. We also compare the data to model calculations of the INS intensities under the same reciprocal space averaging conditions, as shown in Fig.~\ref{dispersion}(f)--(h).

In Figs.~\ref{dispersion}(c) and (d), slices \#1 and \#2 along the $(1-K,2K)$ and $(H,0)$ directions reveal dispersing AE and AO branches emanating from the $\Gamma$ points that highlight their collective character. In slice \#1, the AE and AO dispersions are well-defined up to their respective K-point Dirac crossings at 90 meV and 140 meV. In slice \#2, the AE and AO branches are well-defined up to the M-point with energies of 70 meV and 115 meV. The AE, AO, and Tb dispersions and intensities are consistent with model calculations shown in Fig. \ref{dispersion} (f)--(g).

Traces of broad OE, OO, FE and FO excitations can be seen as high as $\sim$ 250 meV in slices \#1--3, consistent with model predictions of a 230 meV energy cutoff.  However, these modes are  incoherent and overdamped. Incoherent excitations observed in slices \#1-3 form a steep feature centered in the $\Gamma'$ zones which extends from approximately 50 to 180 meV (white rectangle).  These modes are clearly observed in slice \#3, where the $\Gamma$-zone AO and AE modes are suppressed. Slice \#1 also indicates higher energy incoherent modes centered at the K-point and extending from 140--230 meV (yellow rectangle).

In an attempt to capture these features, the current Heisenberg model can be extended to include longer-range interactions within the kagome layer. These interactions can be added with constraints that fix the M-point AO and AE energies while attempting to lower the energy of the OO and OE branches at the $\Gamma'$-point. However, these interactions also introduce unsatisfactory distortions of the magnetic spectra (such as an overall lowering of the magnetic bandwidth) and the broad nature of the high-energy excitations present difficulties in numerical fitting of the extended models. More details of these extended models are provided in Fig.~S1 \cite{SM}.

Rather than pursuing an extended Heisenberg model, a qualitative understanding of the incoherent high-energy excitations can be obtained from constant energy cuts through the excitation spectra, as shown in Fig.~\ref{constE}.  Starting at 60 meV in Fig.~\ref{constE}(a), we see the intense AE and AO conical dispersions that form concentric rings of intensity in the $\Gamma$ zones.  We also see streaks of intensity that extend longitudinally across the $\Gamma'$ zones.  As we move up in energy, the AO branch eventually reaches the zone boundary around 115 meV [Fig.~\ref{constE}(d)] while the longitudinal streak in the $\Gamma'$ zone persists.  At energies above the AO branch cutoff, only the streak through the $\Gamma'$ zones remain [Fig.~\ref{constE}(e)].  As we continue to increase energy, the $\Gamma'$ streaks overlap with scattering intensity that forms at K-points found at the intersection between $\Gamma'$ zones [eg., at ${\bf q}=(\frac{2}{3},\frac{2}{3})$]. Figure~\ref{constE}(h) shows that only these K-point excitations remain at 200 meV.  All magnetic excitations disappear at about 250 meV. For comparison, Heisenberg model calculations of these constant-$E$ cuts can be found in Fig.~S2 \cite{SM}.  We can make several conclusions about these results; (1) the expected spectral weight of the OE/OO/FE/FO branches is associated with incoherent excitations in the $\Gamma'$ zones, (2) the steep incoherent excitation centered at $\Gamma'$ in Fig.~\ref{dispersion} forms the longitudinal streaks in the $\Gamma'$ zones in Fig.~\ref{constE}(c)--(e), and (3) the K-point excitations extend to the highest energies as observed in Fig.~\ref{dispersion}(c).

The broad nature of the high energy modes in reciprocal space suggests that they are localized excitations formed from magnetic clusters within the kagome layer.  We studied both triangular and hexagonal plaquettes (see Fig.~S5 \cite{SM}) and find that the data is accurately described by spin clusters on a single hexagon, as shown in Figs.~\ref{cluster}(a)--(c).  Considering that the time-averaged Mn moments point perpendicular to the kagome layer, the spin patterns shown correspond to the instantaneous transverse components which precess around the $c$-axis. The neutron intensity is estimated by calculating the corresponding static structure factors
\begin{equation}
S({\bf q}) = f^2(q)\Big|\sum_{j=1,6} {\rm e}^{i\phi_j} {\rm e}^{i{\bf q}\cdot {\bf r}_j} \Big|^2
\end{equation}
 as shown in Figs.~\ref{cluster}(d)--(f).  Here, $f(q)$ is the magnetic form factor and $\phi_j$ describes the relative angle of the instantaneous spin direction of the transverse component of spin $j$. Dynamical spin precession averages over the spin direction at each site and only the relative angle between spins around the hexagon is relevant.

The FM cluster in Fig.~\ref{cluster}(a) consists of in-phase precession of the spins ($\phi_{j+1}=\phi_j$), resulting in strong intensity in the $\Gamma$ zones [Fig.~\ref{cluster}(d)]. This corresponds to the AO and AE modes, although these modes are more appropriately described as collective excitations using spin wave theory \cite{Riberolles22}.  The chiral cluster in Fig.~\ref{cluster}(b) is based on the $q=0$ antiferromagnet (AF), which is a chiral magnetic ground state found in kagome systems with nearest-neighbor AF interactions, such as iron jarosite \cite{Inami00, Grohol05}.  In the long-range-ordered $q=0$ AF state, magnetic Bragg peaks appear in $\Gamma'$ zones and are systematically absent in the $\Gamma$ zones.  In the chiral cluster, Mn spins successively rotate by 120$^\circ$ around the hexagon ($\phi_{j+1}=\phi_j + 2\pi/3$) (positive chirality) and result in an intensity pattern consistent with longitudinal streaks observed in the $\Gamma'$ zones [Fig.~\ref{cluster}(e)].  Positive and negative chiralities have the same intensity pattern and are degenerate in the absence of Dzyaloshinskii-Moriya interactions.  Finally, we consider a flat-band cluster, consisting of successive 180$^\circ$ rotations of spins around the hexagon ($\phi_{j+1}=\phi_j+\pi$) shown in Fig.~\ref{cluster}(c). This is the correct Wannier representation of the localized flat band excitations \cite{Bergman08} and generates scattering intensity which is peaked at the K-points surrounded by $\Gamma'$ zones [Fig.~\ref{cluster}(f)].  Thus, the two incoherent modes are associated with chiral and flat band spin correlations on isolated hexagonal plaquettes.  Damped harmonic oscillator analysis of the energy spectra (see Fig.~S3 \cite{SM}) indicate overdamped character with a quality-factor of 1.5--2.

\begin{figure}
\includegraphics[width=1.\linewidth]{cluster_calcs_figure.png}
\caption{\footnotesize Spin patterns for the in-plane component of the Mn moment in (a) ferromagnetic, (b) chiral, and (c) flat-band hexagonal clusters. Numbers label the unique spins in the kagome layer and red circles indicate the direction of spin precession. Panels (d)--(f) show the corresponding static structure factors for the black spins around the hexagonal plaquettes in (a)--(c), respectively.}
\label{cluster}
\end{figure}

Figure \ref{fig:temperature} and Fig.~S4 \cite{SM} compares cuts through excitations in the $\Gamma$ and $\Gamma'$ zones at $T=5$ K and 400 K, just below the Curie temperature of $T_c=420$ K.  Dispersive AO/AE modes and incoherent $\Gamma'$ chiral modes are both observed at 400 K. The intensity of both modes rapidly diminishes above 100 meV as a consequence of softening and mode damping near $T_c$.  The data suggest that the both the $\Gamma$ and $\Gamma'$ responses originate from the Mn sublattice magnetization.

 \begin{figure}
\includegraphics[width=1.0\linewidth]{400K_cuts_compare.png}
\caption{\footnotesize Comparison of constant-energy cuts along the $(-K,2K)$ direction for $T=5$ K (blue) and 400 K (red).  Data are scaled to be proportional to the dynamical susceptibility, $\chi''({\bf q},E)=I({\bf q},E)(1-{\rm exp}(-E/k_{\rm B}T))$.  Cuts are performed with $E_i=250$ meV at energy transfers of 70, 90, 110, and 130 meV (averaged over an energy range of $\Delta E= \pm 5$ meV).  All plots are averaged over reciprocal space ranges of $H=$ [0.9:1.1] and $L=$ [-7:7]. }
\label{fig:temperature}
\end{figure}

\section{Discussion}
The uniaxial FM kagome layers in TbMn$_6$Sn$_6$ display conventional, collective acoustic spin waves when Mn spins precess in-phase within a kagome layer.  However, optical and flat band modes whose precessions are out-of-phase have an incoherent character, remaining localized to a single hexagonal plaquette with heavy damping.  Confinement to a hexagonal plaquette is expected for the flat band excitations of the kagome lattice due to phase cancellation of coherent spin precession on the triangular vertices surrounding a hexagon (see Fig.~S6 \cite{SM}).  Thus, the observed flat bands would, in principle, be consistent with linear spin wave theory were it not for the heavy damping.  Similar to their electronic counterparts, the observation of the flat magnon band in a FM metal is rare: Flat band magnon modes are only reported in the insulating kagome FM Cu[1,3-benzenedicarboxylate] where damping is small \cite{Chisnell15}, but have not previously been reported in FM kagome metals.

The optical modes are not expected to have a localized, chiral character based on the kagome geometry alone and are truly anomalous. Such localization could arise from magnetic frustration. For example, a similar confinement of excitations to hexagonal plaquettes is observed in frustrated pyrochlore antiferromagnets \cite{Lee02}.  Frustration from competing magnetic interactions are possible in FM kagome metals.  For example, itinerant FM order is proposed to coexist with $q=0$ chiral AF order in Co$_3$Sn$_2$S$_2$ \cite{Guguchia20, Lachman20, Zhang21}.  However, FM order is much more robust in $R$Mn$_6$Mn$_6$ and competing intralayer magnetic interactions are not expected to play a significant role.

Clues to this anomalous behavior could come from the itinerant-like character of the optical and flat magnon modes.  For example, the damping of both flat and optical/chiral modes observed here in TbMn$_6$Sn$_6$ and reported in other FM kagome metals, such as YMn$_6$Sn$_6$ \cite{Zhang20} and FeSn \cite{Do21, Xie21}.  Heavy damping is consistent with the coupling between magnons and electrons, providing a window into the correlated electron physics of kagome metals. In TbMn$_6$Sn$_6$, quantum oscillations, photoemission, and DFT calculations indicate massive Dirac fermions approximately 100 meV above $E_F$ \cite{Yin20,Ma21,Li21,Liu21,Gu22} which could provide decay channels for particle-hole excitations at the $\Gamma$-point (intravalley) and K-point (intervalley) (see in Fig.~S8 in the SM).  However, the massive Dirac fermions in TbMn$_6$Sn$_6$ are spin-polarized minority bands with an exchange splitting of $\sim1$ eV \cite{Lee22}, which should strongly suppress any spin-flip scattering channels. Our DFT calculations of the bare electronic spin-flip susceptibility agree with this assessment (Fig~S7 \cite{SM}).  However, the DFT calculations do point to potential sources of Landau damping from trivial electronic bands (see SM \cite{SM}).

In summary, we have observed the magnetic spectral weight associated with optical and flat magnon bands in a FM kagome metal. The excitations are comprised of chiral and flat band spin correlations within a hexagonal plaquette, but are incoherent between plaquettes. Extreme damping of these cluster excitations could be caused by Landau damping from trivial electronic bands, although the localized character of chiral modes may point to more exotic couplings. For example, future work should consider the possibility that staggered flux or orbital currents predicted to circulate on the same plaquettes \cite{Yin22,Mielke22} could impart chirality, localization, and damping to the magnetic excitations.

\section{Acknowledgments} RJM would like to thank Joe Checkelsky for comments. RJM, LK, PPO, BGU, BL, and SXMR's work at the Ames Laboratory is supported by the U.S. Department of Energy (USDOE), Office of Basic Energy Sciences, Division of Materials Sciences and Engineering. TJS, TH, and PC are supported by the Center for the Advancement of Topological Semimetals (CATS), an Energy Frontier Research Center funded by the USDOE Office of Science, Office of Basic Energy Sciences, through the Ames Laboratory.  Ames Laboratory is operated for the USDOE by Iowa State University under Contract No. DE-AC02-07CH11358. A portion of this research used resources at the Spallation Neutron Source, which is a USDOE Office of Science User Facility operated by the Oak Ridge National Laboratory.


text:
\section{Introduction}
Often invoked in the past in case of unexplained peculiarities, the role of magnetic fields has long been a subject of debate in the massive star community. Indeed, on the theoretical side, these stars lack the convective envelopes thought to be responsible for the magnetic dynamo in late-type stars whereas, on the observational side, line dilution due to bright companions, spectral contamination by emission, and the relative rarity and breadth of spectral lines limit the detectability of magnetic Zeeman signatures in hot stars. At the beginning of this century, only one half-dozen hot stars (i.e. with spectral types of at least B2) had been confidently detected to be magnetic; the earliest of this sample was only of spectral type B1.5. In the last decade, however, the situation has improved dramatically thanks to the advent of efficient spectropolarimeters attached to mid- and large-size telescopes (such as e.g. FORS at the ESO VLT or ESPaDOnS at the CFHT). These instruments led to the detection of magnetic fields in about 30 hot stars (for a few recent cases: \citealt{riv10,hub11,pet11,gru12}).

The role of magnetic fields in massive stars has thus been re-assessed. While \citet{mae03} focused on understanding their effect on the structure and evolution of massive stars, \citet{bab97a} had unveiled their impact on the stellar winds and therefore on the X-ray emission of massive stars. In a seminal paper, they proposed a semi-analytic model for magnetically channeled wind shocks wherein a strong dipole magnetic field guides the stellar winds from two opposite hemispheres towards the magnetic equator where the flows collide: the gas is heated up to high temperatures, leading to hard X-ray emission. Such channelling can naturally lead to a shock-heated equatorial region.

This model was first developed and applied to explain the peculiar X-ray emission of the magnetic and chemically peculiar Bp star IQ\,Aur. Following the work of \citet{sta96}, it was extended to explain the peculiarities of the O7 star $\theta^1$\,Ori\,C \citep{bab97b}. Self-consistent numerical magnetohydrodynamic models by \citet{udd02} further explored the wind-field interaction in hot stars, ultimately establishing $\theta^1$\,Ori\,C as a prototype for which the 15d modulations of the \ha, UV, and X-ray emissions are well explained by a magnetic oblique rotator model. Furthermore, \citet{udd06} have shown that in the presence of significant stellar rotation, centrifugal support can lead to a stable accumulation of material near the magnetic equator. The best example here is the Bp star $\sigma$\,Ori\,E, for which the magnetic energy is so strong compared to the stellar wind energy that the magnetosphere extends to many stellar radii \citep{tow05}. For both cases, the X-ray emission appears very peculiar, as it is altogether bright and hard. Indeed, ``normal'' O-type stars have relatively soft X-ray emission ($kT\sim0.2-0.7$\,keV) of a well-defined strength (the so-called canonical relation, \loglxlb $\sim10^{-7}$, see \citealt{naz11,naz09} and references therein). Both $\theta^1$\,Ori\,C and $\sigma$\,Ori\,E therefore seem to comply reasonably well with model expectations for magnetically confined winds.

It must however be noted that, while bright and hard X-ray emission may be considered as a good hint for the presence of magnetic fields, the contrary is not necessarily true. Indeed, some magnetic objects display soft X-ray emission: this is notably the case of Of?p stars (see e.g. \citealt{naz10,naz12}) and some B stars (see \citealt{gud09,osk11} and references therein). While the overall X-ray luminosity of these objects follows rather well the expected values, it is clearly not the case of their plasma temperature \citep{osk11}, a fact which currently remains unexplained.

The Chandra X-ray observatory has recently performed a large-scale survey of the Carina nebula \citep[CCCP, see ][]{tow11}. This led to the detection of $>$14\,000 X-ray sources \citep{bro11}, among which a majority of PMS objects and $>$100 massive, hot stars \citep{naz11}. Amongst the latter ones, we selected those that display unusually hard and bright X-ray emission, hence were good candidates for harbouring magnetically-confined winds. Using FORS, we perform a small spectropolarimetric survey of these objects, with the aim of finding new cases of magnetically-confined winds. This paper reports on the results of this campaign. It first presents the observations (Sect. 2), then the results (Sect. 3), to finally summarize and conclude in Sect. 4.

\section{Observations}

\subsection{Observational Strategy}

A total of eight fields in the Carina nebula were observed in spectropolarimetric mode with the FORS2 instrument \citep{app98} of the ESO VLT during a visitor run in March 2011. The choice of the main targets was initially made on the basis of their X-ray emission \citep[CCCP, ][and references therein]{naz11}. We selected bright\footnote{$V<11.$\,mag} early-type stars with high X-ray luminosity (Coll228-68, HD\,93250, Tr16-5, 11, 22 and 64), with a strong hard X-ray component (HD\,93250, HD\,93501, LS\,1865, Coll228-68, Tr16-5, 10, 11, 22 and 64) or with known long-term X-ray variability (HD\,93250, Tr16-22). Those characteristics are compatible with those expected from magnetically-confined winds, explaining our choice for this investigation.

Four fields studied only a single target, because of the lack of nearby massive stars with similar magnitudes: HD\,93250, HD\,93501, HD\,93190, and Tr16-22. The four other fields were taken in multi-object mode, since such adequate neighbours existed: the first one contains Coll228-66, Coll228-68, and HD\,93097, the second one Tr16-10, 11 and 14, the third one LS\,1865, LS\,1853, HD\,305524, and HD\,305534, and the last one Tr16-2, 4, 5, 13, 15, 64, and 115.  The 11 additional objects are field stars which do not show particularly exceptional X-ray characteristics.

In total, 21 OB stars were thus investigated. Our observing run was allocated 10\,hours, that were distributed over three consecutive nights (2011  March 12--14). This allowed us to observe about half of our targets on two different nights.

During our observing run, the blue CCD (a mosaic composed of two 2k$\times$4k E2V chips) read without any binning, a slit of 1\arcsec\,\footnote{The second night suffered from a bad seeing, which forced us to lengthen the exposure durations, and even use a 1.2\arcsec slit for one of the fields (that centered on Tr16-22). A consequence of this was a change in extraction radius during the reduction (see below).} and the 1200B grating ($R\sim 1400$) were used. We adopted for each field an observing sequence of 8 or 16 subexposures with retarder waveplate positions of $+45^{\circ}$, $+45^{\circ}$, $-45^{\circ}$, $-45^{\circ}$, $+45^{\circ}$, $+45^{\circ}$, $-45^{\circ}$, $-45^{\circ}$,...  Duration of individual subexposures ranges from 10 to 300s, depending on the targets' luminosity and the variable seeing. It was therefore adjusted several times during the run, to remain close, but below, the saturation level (spectropolarimetric investigations requiring the highest possible signal-to-noise ratio).

\subsection{Data Reduction}

In recent years, there has been some debate in the literature about the reality of some magnetic field detections claimed using the FORS instruments. To settle the debate, \citet{bag12} has performed an overall reduction of all spectropolarimetric FORS datasets. They notably show the implications of some data reduction choices and we took into account the considerations made in that work.

Using the observing sequence described above, the normalized Stokes $V/I$ profile can be determined, as well as a diagnostic ``null''  profile \citep{don97,bag09}:\\
\begin{equation}
\begin{array}{rcl}
\pv = V/I &=& {1 \over 2 N} \sum\limits_{j=1}^N \left[
\left(\frac{\fo - \fe}{\fo + \fe}\right)_{-45^\circ} -
\left(\frac{\fo - \fe}{\fo + \fe}\right)_{+45^\circ}\right] \\[2mm]
\nv = N/I &=& {1 \over 2 N} \sum\limits_{j=1}^N (-1)^{(j-1)}\left[
\left(\frac{\fo - \fe}{\fo + \fe}\right)_{-45^\circ} -
\left(\frac{\fo - \fe}{\fo + \fe}\right)_{+45^\circ}\right]\; ,\\
\end{array}
\label{Eq_V_and_N}
\end{equation}
where \fo\ and \fe\ are the counts recorded for the parallel and perpendicular beams, respectively, $I$ is the simple sum of all \fo\ and \fe\ spectra, and $N$ is the number of pairs of subexposures (i.e. sets of $\{-45^\circ$, $+45^\circ \}$ observations). To get the \fo\ and \fe\ for each target, the data were corrected for bias but not flat-fielded, since there is no evidence that flatfielding improves the magnetic field determination \citep{bag12}. Furthermore, the spectra were extracted using a simple aperture extraction since,  in practice, optimal extraction applied to high signal-to-noise ratio leads to noisier spectra \citep[e.g. ][]{bag12}. The aperture extraction radius was fixed to 12\,px for data taken on the first and third night, and between 15 and 24\,px (depending on the actual seeing) for data taken during the second night. Note that sky background was subtracted for each spectrum\footnote{The background subtraction in IRAF (see below) used an average of small sky regions surrounding the star (usually $-35:-25$ and $+25:+35$, or $-55:-45$ and $+45:+55$ when the seeing was bad or a companion existed - only one of these regions being used if the star was too close to a slit edge).}. Finally, wavelength calibration was made by considering arc lamp data taken at only one retarder waveplate position (in our case, $-45^\circ$) for all science frames of a given field, whatever the position of the retarder waveplate actually is for them. As shown by \citet{bag09}, this minimizes spurious polarization signals. Note that, indeed, wavelength calibration remains specific to the actual position of the spectra on the CCD, i.e. it is different from one object to the next as well as for the parallel and perpendicular beams. In all cases, the final, calibrated spectra have 5813\,px with a step of 0.43\AA\ (close to the dispersion value) and a start wavelength of 3500.215\AA, as recommended for this grating.

Whenever it was possible, we reduced our data using two different procedures. For the data obtained in ``fast mode" (i.e. single target), we used both the FORS pipeline\footnote{ftp://ftp.eso.org/pub/dfs/pipelines/fors/fors-pipeline-manual-4.1.pdf} and standard IRAF routines\footnote{$imcombine$ (for master bias creation), $apall$ (for spectrum extraction), and $identify$ (for wavelength calibration, using the wavelength list from the FORS manual which is available at http://www.eso.org/sci/facilities/paranal/instruments/fors/doc/).}. However, the FORS pipeline only accepts datasets with a pair of frames obtained with the retarder waveplate at $+45^{\circ}$ and $-45^{\circ}$, or a quadruplet with retarder waveplate at the position angles $+45^{\circ}$, $-45^{\circ}$, $+135^{\circ}$, and $+225^{\circ}$. Since we adopted a different observing strategy (8 or 16 subexposures taken at two angles, rather than a single pair, see above), we had to combine our 8 or 16 subexposures into 4 or 8 pairs, to reduce these pairs separately, and finally to make use of external procedures to combine the pipeline products into a single result for each 8 or 16 subexposure dataset. Moreover, for polarimetric data obtained in multi-object mode, we could use only IRAF routines, because the FORS pipeline fails to associate the pairs of beams pertaining to individual objects. In the end, the spectra reduced within IRAF appeared very similar to those reduced with the FORS pipeline, when they existed, though the former ones seem slightly less noisy. This is explained by bad seeing, and by the position of the spectra, which, even in ``fast'' mode, were found very close to gaps between slitlets. These problems could be better dealt with through manual fine tuning of extraction and background substraction, than with the pipeline semi-automatic mode. Nevertheless, we decided to use also the FORS pipeline to check our results in two independent ways.

In a final step, the longitudinal magnetic field was estimated by minimising the expression \citep{bag02}:
\begin{equation}
\chi^2 = \sum_i \frac{(y_i - \bz\,x_i - a)^2}{\sigma^2_i}
\label{Eq_ChiSquare}
\end{equation}
where $i$ takes a different value for each spectral point, $a$ is a constant accounting for possible spurious residual polarization in the continuum, $y_i$ is either \pv\ or \nv\ at the wavelength $\lambda_i$, and $x_i = -g_\mathrm{eff}\ \cz \ \lambda^2_i\ 1/I_i\ (\mathrm{d}I/\mathrm{d}\lambda)_i$, with $\cz \simeq 4.67 \times 10^{-13}$\,\AA$^{-1}$\,G$^{-1}$ and $g_\mathrm{eff}$ the effective Land\'{e} factor, assumed to be 1 near the Balmer lines of hydrogen and 1.2 elsewhere.

In this equation, $(\mathrm{d}I/\mathrm{d}\lambda)_i$ was evaluated through $(I_{i+1}-I_{i-1})/(\lambda_{i+1}-\lambda_{i-1})$, as by \citet{bag02}, while the error on \pv\ is estimated through $\sigma_i=\frac{1}{2\sqrt{N}} \, \frac{\sigma(f)_i}{f_i}$ \citep[see eq. A6 of][]{bag09}. Following the FORS pipeline manual\footnote{ftp://ftp.eso.org/pub/dfs/pipelines/fors/fors-pipeline-manual-4.1.pdf}, the error on each spectrum (\fo\ or \fe) amounts to $\sigma(f)_i=\sqrt{f_i/g}$ where $f_i$ is the flux at the wavelength $\lambda_i$ and $g$ is the gain in electrons per ADU (0.53 or 0.58, depending on the CCD chip) - the read-out noise is neglected because of its small value compared to the high level of the stellar flux. Note that the FORS pipeline accounts for the error introduced by background subtraction and by the extraction, whereas for data reduced with IRAF routines we considered only the contribution due to photon-noise.

A few additional remarks must be made at this stage. First, we have estimated the magnetic field for the ``raw'' \pv\ and \nv\ profiles, and for the same \pv\ and \nv\ profiles rectified using a Fourier filter of width 300px \citep{bag12}. That correction results only in a minimal difference (about 15G on average), well within the error bars (typically 80G), confirming the results of \citet[see their Table 2]{bag12}. Second, some values of the \pv\ and \nv\ profiles were discarded: the 3\% bluest and 3\% reddest of the wavelength range with recorded data, as well as the spectral bins where \nv\ deviates more than 3$\sigma$ from zero (and the two adjacent bins), as in \citet{bag12}. In addition, we did not consider the spectral bins with $x_i$ outside the range $\{-1\times 10^{-6}$, $+1\times 10^{-6} \}$, as they were small in number and erratic in value but dominating the slope determination. Third, we derived magnetic field values for the full spectral range and for selected spectral windows (avoiding continuum regions, or regions affected by ISM/nebular features, see online material in Appendix). Again, this results in small changes (about 45G on average), within the error bars. Finally, \citet[and references therein]{bag12} found several hints that the error bars of the magnetic field determinations made using the FORS instrument were underestimated. Therefore, the error on \bz, derived as is usual from the linear interpolation, was corrected following the method of \citet{bag12}, to take into account the actual scattering around the fitted line. This correcting factor is $\sqrt{\chi^2_{\rm min}/\nu}$, where $\chi^2_{\rm min}$ is the minimum value of the $\chi^2$ derived as in Eq.~(\ref{Eq_ChiSquare}) and $\nu$ is the number of degrees of freedom of the system (i.e., the number of spectral points minus two).

\section{Results}

The derived values of the longitudinal fields for all targets are listed in Table \ref{result}, which also gives for comparison the \nz\ values found from the null profiles. For comparison purposes, values with or without rectification, for the full spectral range or considering only specific spectral windows are given. Spectral types and X-ray luminosity, reproduced from \citet{naz11}, are also shown. Several conclusions can now be drawn.

First, even though the FORS pipeline results are about 10-15\% noisier then those obtained with manual data reduction, all results are consistent with each other (i.e. the conclusion of a detection or a non-detection is found with both methods), and the field values are within 2$\sigma$ error bars of each other. Second, the diagnostic \bz\ derived from the null profile are consistent with zero (at the 3$\sigma$ level in the worst cases, most often at the 1$\sigma$ level), as is expected, confirming the quality of our data reduction. Finally, and most importantly, two stars shows a magnetic field detection at the 6$\sigma$ level and with a strength $>$300\,G, two features suggesting a true detection \citep{bag12}. The other 19 stars have \bz\ compatible with zero at the 3$\sigma$ level, as do the diagnostic profiles.

\begin{figure*}
\includegraphics[width=8cm]{stokesv_22n2_rect.ps}
\includegraphics[width=8cm]{stokesv_13_rect.ps}
\caption{Normalized Stokes V (in selected spectral windows, with rectification) as a function of x ($= -g_\mathrm{eff}\ \cz \ \lambda^2\ 1/I\ (\mathrm{d}I/\mathrm{d}\lambda)$, see Eq.~(\ref{Eq_ChiSquare})) for the two positive cases of Tr16-22 (left) and Tr16-13 (right). In these panels, the slope is directly linked to the intensity of the longitudinal field. The best-fit lines, derived from the least-squares analysis, are drawn with thick dashed lines, for both the actual Stokes V parameter (top panels) and the diagnostic null profile (bottom panels). In the latter case, or more generally in the case of non-detections, these lines should be horizontal.  The line thickness is larger than the 1$\sigma$ error bar on the constant term, while the two thin dotted lines show lines with slopes at $\pm 1 \sigma$ from the best-fit value.}
\label{stokesv}
\end{figure*}

\begin{table*}
\begin{center}
\caption{Field measurements for all stars of our sample, for both Stokes $V$ profiles and their related diagnostic null profiles. For fields with only one target, the results found from the FORS pipeline are also quoted at the bottom of the table. A star after the name indicates the main targets, those displaying X-ray peculiarities, and boldfaces shows the 6$\sigma$ detections.}
\label{result}
\begin{tabular}{l l c | r r | r r | r r }
            \hline\hline
Target & Sp. Type & \loglxlb & \multicolumn{2}{|r|}{Full range, non-rect.} & \multicolumn{2}{r|}{Spectral windows, non-rect.} & \multicolumn{2}{r}{Spectral windows, rectified}\\
 & & & \bz & \nz & \bz & \nz & \bz & \nz \\
 & & & (G) & (G) & (G) & (G) & (G) & (G) \\
\hline
\multicolumn{9}{c}{IRAF} \\
Coll228-66           & O9.5V   & $-7.5$  & $-88\pm43$ & $14\pm42$  & $-72\pm45$ & $-6\pm44$  & $-63\pm45$  & $-6\pm44$  \\
Coll228-68$\star$    & B1Vn    & $-6.5$  & $31\pm56$  & $-31\pm51$ & $60\pm58$  & $-48\pm55$ & $78\pm58$   & $-46\pm55$ \\
HD\,93097            & B0Vn    & $<-7.9$ & $-57\pm44$ & $-3\pm44$  & $-89\pm47$ & $31\pm48$  & $-80\pm46$  & $1\pm47$   \\
HD\,93190$\star$     & B0IV:ep & $-7.8$  &$-120\pm93$ & $-10\pm90$ &$-154\pm120$& $-25\pm120$& $-199\pm128$& $26\pm125$ \\
HD\,93250-n1$\star$  &O4III(fc)& $-6.4$  & $96\pm78$  &$123\pm75$  & $90\pm85$  & $94\pm83$  & $58\pm82$   & $96\pm82$  \\
HD\,93250-n2$\star$  &         &         & $25\pm69$  & $49\pm69$  & $17\pm74$  & $71\pm72$  &  $8\pm73$   & $43\pm71$  \\
HD\,93501-n1$\star$  & B1.5III:& $-6.9$  &  $9\pm74$  & $201\pm72$ & $-3\pm81$  & $267\pm79$ & $-3\pm81$   & $272\pm78$ \\
HD\,93501-n2$\star$  &         &         & $60\pm70$  & $-4\pm69$  & $12\pm76$  & $-58\pm75$ & $18\pm76$   & $-66\pm74$ \\
HD\,305524-n1        & O7V((f))& $-7.3$  & $-50\pm71$ & $-96\pm69$ & $-65\pm79$ & $-72\pm75$ & $-16\pm76$  & $-79\pm75$ \\
HD\,305524-n2        &         &         & $-120\pm57$& $-13\pm51$ & $-97\pm68$ & $-3\pm58$  & $-55\pm64$  & $1\pm58$   \\
HD\,305534-n1        & B0.5V+B1V:& $-7.9$& $-4\pm78$  & $-124\pm78$& $14\pm87$  & $-197\pm87$& $37\pm86$   & $-199\pm87$\\
HD\,305534-n2        &           &       & $-19\pm59$ & $-39\pm56$ & $-46\pm66$ & $-86\pm65$ & $-63\pm65$  & $-88\pm65$ \\
LS\,1853-n1          & B1Ib    & $<-9.0$ & $41\pm60$  & $12\pm60$  & $48\pm64$  & $-6\pm62$  & $48\pm63$   & $-9\pm62$  \\
LS\,1853-n2          &         &         & $-6\pm52$  & $-103\pm51$& $18\pm55$  & $-116\pm55$& $24\pm54$   & $-118\pm55$\\
LS\,1865-n1$\star$   &O8.5V((f))& $-6.9$ & $-72\pm67$ & $-83\pm67$ & $-41\pm72$ & $-84\pm69$ & $-56\pm72$  & $-86\pm70$ \\
LS\,1865-n2$\star$   &          &        & $-70\pm52$ & $64\pm51$  & $-98\pm56$ & $92\pm52$  & $-112\pm56$ & $95\pm52$  \\
Tr16-2               & B1V     & $<-8.0$ & $120\pm62$ & $33\pm61$  & $82\pm66$  & $62\pm65$  & $55\pm65$   & $51\pm65$  \\
Tr16-4               & B1V     & $<-7.9$ & $1\pm59$   & $113\pm60$ & $29\pm63$  & $100\pm62$ & $21\pm63$   & $101\pm62$ \\
Tr16-5$\star$        & B1V     & $-6.4$  & $86\pm76$  & $-115\pm74$& $42\pm85$  & $-131\pm83$& $19\pm85$   & $-128\pm83$\\
Tr16-10$\star$       & B0V     & $-7.0$  & $124\pm102$& $60\pm100$ & $121\pm108$& $116\pm108$& $153\pm108$ & $114\pm108$\\
Tr16-11$\star$       & B1.5V   & $-6.4$  &$-104\pm201$&$-177\pm205$&$-155\pm225$& $168\pm230$& $-145\pm229$& $175\pm230$\\
Tr16-13              & B1V     & $-7.8$  & $\bld{490\pm70}$ & $-7\pm64$  & $\bld{456\pm76}$ & $13\pm71$  & $\bld{393\pm73}$  & $16\pm71$  \\
Tr16-14              & B0.5V   & $-8.0$  &$-143\pm166$& $161\pm165$& $91\pm186$ & $433\pm191$& $120\pm187$ & $418\pm190$\\
Tr16-15              & B0V     & $-7.9$  & $18\pm68$  & $24\pm68$  & $-1\pm73$  & $65\pm73$  & $-6\pm73$   & $67\pm72$  \\
Tr16-22-n1$\star$    & O8.5V   & $-6.3$  & $\bld{-386\pm70}$& $3\pm68$   &$\bld{-433\pm73}$ & $54\pm73$  & $\bld{-454\pm72}$ & $45\pm74$  \\
Tr16-22-n2$\star$    &         &         & $\bld{-450\pm81}$& $-28\pm76$ &$\bld{-578\pm88}$ & $-94\pm83$ & $\bld{-604\pm87}$ & $-88\pm82$ \\
Tr16-64$\star$       & B1.5Vb  & $-6.1$  & $-43\pm55$ & $17\pm53$  & $-49\pm58$ & $41\pm57$  & $-65\pm57$  & $43\pm57$  \\
Tr16-115             & O9.5V   & $-7.2$  & $-15\pm67$ & $38\pm66$  & $-31\pm69$ & $23\pm67$  & $-38\pm68$  & $22\pm67$  \\
\hline
\multicolumn{9}{c}{FORS pipeline} \\
HD\,93190$\star$     & B0IV:ep & $-7.8$  & $-272\pm110$ & $-23\pm104$ & $-192\pm141$&$-45\pm138$& $-154\pm140$& $2\pm140$  \\
HD\,93250-n1$\star$  &O4III(fc)& $-6.4$  & $102\pm91$   & $120\pm87$  & $95\pm98$   & $92\pm96$ & $68\pm95$   & $78\pm95$  \\
HD\,93250-n2$\star$  &         &         & $-128\pm137$ & $82\pm116$  & $-250\pm144$& $93\pm126$& $-309\pm139$& $47\pm122$ \\
HD\,93501-n1$\star$  & B1.5III:& $-6.9$  & $13\pm85$    & $200\pm81$  & $-3\pm93$   & $295\pm90$& $22\pm92$   & $309\pm89$ \\
HD\,93501-n2$\star$  &         &         & $90\pm84$    & $41\pm78$   & $77\pm93$   & $3\pm87$  & $84\pm91$   & $-1\pm86$  \\
Tr16-22-n1$\star$    & O8.5V   & $-6.3$  & $\bld{-463\pm83}$  & $-43\pm76$  & $\bld{-575\pm91}$ & $3\pm81$  & $\bld{-620\pm89}$ & $-5\pm82$  \\
Tr16-22-n2$\star$    &         &         & $\bld{-483\pm92}$  & $29\pm86$   & $\bld{-604\pm101}$& $0\pm97$  & $\bld{-632\pm99}$ & $-8\pm98$  \\
\hline
\end{tabular}
\end{center}
{\footnotesize {1. Comparison of names in Simbad or \citet{naz11} and those of this paper: CPD -59255B = Coll228-66, Cl* Trumpler 14 MJ 449 = LS\,1865, Cl* Trumpler 14 MJ 366 = LS\,1853, Cl* Trumpler 16 MJ 506 = Tr16-2, Cl* Trumpler 16 MJ 466 = Tr16-4, Cl* Trumpler 16 MJ 427 = Tr16-5, Cl* Trumpler 16 MJ 327 = Tr16-10, Cl* Trumpler 16 MJ 289 = Tr16-11, Cl* Trumpler 16 MJ 339 = Tr16-13, Cl* Trumpler 16 MJ 357 = Tr16-14, Cl* Trumpler 16 MJ 372 = Tr16-15, Cl* Trumpler 16 MJ 496 = Tr16-22, Cl* Trumpler 16 MJ 554 = Tr16-115.\\
2. Note that the typical \loglxlb\ in Carina Chandra data is $-7.2$ \citep{naz11}.}}
\end{table*}

\subsection{Tr16-22}
With the Chandra and XMM observations of the Carina nebula, Tr16-22 has been singled out as a peculiar object. Indeed, it is strongly overluminous in X-rays \citep[and references therein]{eva04,ant08,naz11}. Moreover, this high-energy emission is unusually hard for an O-type star, with a component at 2\,keV of a strength similar to that of the soft component at 0.7\,keV \citep{naz11}. Furthermore, the X-ray luminosity varies by a factor of 30 \citep{com11}, whereas the X-ray emission of O-stars is known to be rather stable.

\begin{figure}
\includegraphics[width=8cm]{feros22_new.ps}
\caption{Archival FEROS spectrum of Tr16-22, showing the narrow metallic lines, the P cygni profile of He\,{\sc i}\,5876\AA, and the strong \ha\ emission. Note that the narrowest emission component in \ha\ is nebular in origin (other nebular lines, e.g. those of [N\,{\sc ii}], are indeed detected). }
\label{feros}
\end{figure}

Radial velocity variations have been reported for this object from low-resolution monitorings \citep{com11,wil11}, suggesting the star to be a binary. With the information at hand before this study, X-rays linked to a wind-wind collision might at first seem a good explanation for the observed high-energy characteristics. However, the observed large overluminosity ($>$1\,dex) is neither observed nor theoretically expected in late-type O+OB colliding-wind binaries \citep[for examples of O+O binary X-ray observations, see][]{naz09,naz11}. Moreover, the reported spectral changes, observed with low signal-to-noise \citep{com11}, could equally well be compatible with line profile variations. In view of these facts and our spectropolarimetric result (Table \ref{result} and Fig. \ref{stokesv}), magnetically-confined winds now appear as a much better explanation \citep{eva04,naz11}. A monitoring is now needed to further constrain the stellar, wind, and magnetospheric properties which would enable us to perform a full modelling of this star. In this context, we have examined the single archival FEROS spectrum of Tr16-22. This spectrum appears quite similar to that of the sharp-lined magnetic O9.5IV star HD\,57682 \citep{gru09}, though the lines of Tr16-22 seem even narrower than those of HD\,57682. A long period for Tr16-22 would not be surprising, as an intense field is able to brake the rotation of massive stars with strong stellar winds \citep{udd09}. Indeed, many magnetic O stars are slow rotators (see e.g. the extreme cases of 538d for HD\,191612 or even 55yrs for HD\,108, \citealt{naz10}). The FEROS spectrum also reveals a strong emission component in \ha\ and a P Cygni profile for He\,{\sc i}\,5876\AA\ (Fig. \ref{feros}), two features comparable to those of the strongly confined $\theta^1$\,Ori\,C and of the Of?p stars. This reinforces our suspicion that Tr16-22 shows several evidence for the presence of magnetically confined winds.

\subsection{Tr16-13}
This object displays a clear magnetic detection (Table \ref{result} and Fig. \ref{stokesv}), but it could only be observed once, so that we have no independent confirmation of our result. Nevertheless, the obtained field values clearly differ between the actual Stokes V profile and the null profile, leaving little doubt on the detection. Spectra of targets observed at the edge of the field of view may be slightly depolarized \citep{bag02}, but this problem would decrease the absolute field value, rather than causing a spurious detection, thus we are confident in our result.

Tr16-13 is not known to exhibit any peculiarity (e.g. spectroscopic, photometric, or radial velocity variations), but the literature concerning this star is very sparse (e.g. no high-resolution spectra appear to exist). The FORS spectrum does not show pronounced He chemical peculiarities, but this is not unusual amongst early B-type magnetic stars (see e.g. HD\,66665, \citealt{pet11}). The spectral lines of Tr16-13 appear somewhat broad (FWHM of about 170\,\kms), but at the limit of the FORS resolution: convolved to the same resolution, the spectrum of $\tau$\,Sco, a sharp-lined, magnetic B0.5V star, appears similar but with narrower lines. However, contrary to that object, Tr16-13 emits few X-rays. Line profile variations may be expected for this star, and follow-up investigations are thus clearly needed to pinpoint its physical properties.

\subsection{HD\,93250}
The non-detection of a magnetic field may naively appear particularly puzzling in the case of HD\,93250. This object indeed displays non-thermal radio emission \citep{lei95}. Such radiation is known to arise from relativistic electrons accelerated at the shocks between the two winds in a binary system \citep{van}. Moreover, HD\,93250 displays large X-ray variations \citep{rau09}, which reinforces the suspicion of a binary nature. Though radial velocity variations have never been observed up to now \citep[see][and references therein]{rau09,wil11}, a close companion was recently detected using interferometry \citep{san11}. HD\,93250 now appears composed of two similar stars, whose winds collide, in a long-period, eccentric orbit. A magnetic field is indeed required in this system, to explain the synchrotron radio emission of the system. However, the nature and origin of the field - whether associated directly to the star(s) or produced within the shock - is unclear. Considering the first case, it is true that our observations do not detect any organized field at the photosphere, but it must be noted that a low inclination of the star's magnetic axis or its binary nature could impact on our detection capabilities. Moreover, the derived 3$\sigma$ limit on the longitudinal field of HD\,93250 is about 300G, and no theoretical expectation has been published to which we could compare this value. Depending on the exact stellar and wind properties (which are not known with precision), as well as local characteristics of the wind collision zone (equally unknown), which may produce a field amplification \citep{bel04,gon12}, it may well be that our non-detection of an organized magnetic field is not at odds with the presence of the $\sim$1\,G field needed at the collision zone to generate some non-thermal radio emission.

\section{Summary and Conclusion}
Following the detection of hard and bright X-rays in a handful of early-type stars in the Carina nebula, we performed a small spectropolarimetric survey of 21 massive stars in the region using the FORS instrument. Previously, only 7 OB stars in Carina (HD\,93128, HD\,93129B, HD\,93204, HD\,93403, HD\,93843, HD\,303308, HD\,303311) had been examined by \citet{hub11}, alas with no clear detection and only one 3$\sigma$ detection \citep[which is a low detection level for FORS data, see ][]{bag12}.

With this paper, we have quadrupled the number of Carina massive stars searched for magnetic fields. For 19 objects, we could put limits of a few hundreds gauss on magnetic fields, but two objects led to a clear (6$\sigma$) detection: Tr16-22 ($\bz \sim 500$\,G) and Tr16-13 ($\bz \sim 450$\,G).

Though the late-O star Tr16-22 could have a companion, its hard, variable X-ray emission appears too bright to originate in a wind-wind collision inside the binary. The magnetic field detection (confirmed by two independent observations) rather supports a scenario where magnetically-confined winds collide at the magnetic equator, thereby producing copious hard X-ray emission. An independent confirmation comes from its strong \ha\ emission, similar to that observed in other magnetic cases (the Of?p stars and $\theta^1$\,Ori\,C). With its peculiar X-ray properties, Tr16-22 may thus be seen as an analog to the famous $\theta^1$\,Ori\,C, the best example amongst magnetic O-stars for a good match with theoretical expectations, especially in the X-ray range \citep{gag05}: as a potential second case, Tr16-22 appears as a key in our understanding of confined winds. A full monitoring of this star (X-ray, spectroscopy, spectropolarimetry,...) is thus urgently needed.

Tr16-13 is a poorly known star with no known peculiarities and rather faint X-ray emission. Certainly, more investigation is now needed to ascertain the impact of magnetic fields on its evolution and on its stellar wind. Again, a dedicated monitoring is required to determine the physical parameters of the star as well as the magnetospheric properties, needed before any modelling can take place.

\section*{Acknowledgments}
The authors acknowledge the referee, I.D. Howarth, for his useful comments and promptness to review the paper. We also thank A. David-Uraz, C. Folsom, M. Shultz for useful comments and contributions for a follow-up ESO proposal.
YN acknowledges support from the Fonds National de la Recherche Scientifique (Belgium), the Communaut\'e Fran\c caise de Belgique, the PRODEX XMM and Integral contracts, and the `Action de Recherche Concert\'ee' (CFWB-Acad\'emie Wallonie Europe). VP acknowledges support from the Fonds Qu\'ebecois de la Recherche sur la Nature et les Technologies. GAW ackowledges Discovery grant support from the Natural Science and Engineering Research Council of Canada. ADS and CDS were used for preparing this document.