# Thermal decay of two-spinon bound states in quasi-2D triangular
antiferromagnets

I. L. Pomponio1, E. A. Ghioldi2, C. J. Gazza1, L. O. Manuel1, and A. E.
Trumper1 1 Instituto de Física Rosario (CONICET) and Facultad de Ciencias
Exactas, Ingeniería y Agrimensura, Universidad Nacional de Rosario, 2000
Rosario, Argentina  2 Department of Physics and Astronomy, University of
Tennessee, Knoxville, Tennessee 37996-1200, USA

###### Abstract

We analyze the temperature evolution of the anomalous magnetic spectrum of the
spin-1/2 triangular quantum Heisenberg antiferromagnet, which is proximate to
a quantum critical point. Recently, its low energy excitations have been
identified with two-spinon bound states, well defined in an ample region of
the Brillouin zone. In this work, we compute the thermal magnetic spectrum
within a Schwinger boson approach, incorporating Gaussian fluctuations around
the saddle-point approximation. In order to account for a finite Néel
temperature $T_{N}$, we incorporate an exchange interaction between triangular
layers. As temperature rises, the dispersion relation of the two-spinon bound
states, representing single-magnon excitations, remains unchanged but becomes
mixed with the thermally activated spinon continuum. Consequently, a crossover
occurs at a temperature $T^{*}\simeq 0.75T_{N}$, defining a terminated
Goldstone regime between $T^{*}$ and $T_{N}$, where only the Goldstone modes
survives as well-defined magnon excitations, up to the Néel temperature. Our
results support the idea that the fractionalization of magnons near a quantum
critical point can be extended to more realistic quasi-2D frustrated
antiferromagnets.

##  I Introduction

The understanding of the organizing principles [1] governing frustrated
antiferromagnets (AF) remains a debated topic in condensed matter physics [2].
The pursuit of quantum spin liquids —highly entangled quantum states without
classical analogs and characterized by fractional spinon excitations— has been
a central focus since the proposal of the resonant valence bond state for the
triangular antiferromagnet [3, 4, 5, 6, 7, 8, 9]. This quest has driven both
theoretical and experimental efforts, fueled by insights from high-temperature
superconductors [10] and concepts from the fractional quantum Hall effect
[11]. Advancements in numerical and analytical techniques, applied to a
variety of frustrated models, have played a relevant role in validating zero-
temperature quantum phase diagrams along with their corresponding low-lying
magnetic excitations [12]. Additionally, the synthesis of new compounds [13]
and the refinement of experimental techniques, such as inelastic neutron
scattering (INS) [14], have confirmed the experimental realization of
predicted exotic magnetic states of matter. In the current era of quantum
materials [15, 16], it is crucial to integrate this knowledge to enable
theoretical predictions of experiments under diverse conditions in an even
more controlled manner.

In recent years, the synthesis of triangular antiferromagnets has notably
increased, giving rise to numerous candidates for spin liquids [17, 18].
However, only few of them are perfectly equilateral, disordered-free, or
faithfully realize a Heisenberg-like Hamiltonian [19]. Instead of directly
pursuing highly exotic magnetic states solely characterized by the absence of
specific features (such as no ordering down to the lowest temperature or
absence of magnon excitations), a more productive strategy has emerged.
Specifically, the study of magnetically ordered or disordered compounds near a
quantum critical point has proven to be fruitful. A notable example is the
well-studied compound Ba3CoSb2O9 that realizes an effective spin-1/2
triangular antiferromagnet with XXZ model interactions, exhibiting a 120° Néel
order below $T_{N}$ = 3.8K [20, 21, 22, 23, 24]. In INS experiments conducted
below $T_{N}$, an unusual spectrum was observed, comprising low-lying energy
collective excitations (magnonic-like) coexisting with two dispersive continua
at higher energies across the entire Brillouin zone (BZ). The semiclassical
theory (large-$S$) [23, 25, 26] falls short in explaining this spectrum,
leading to speculation that Ba3CoSb2O9 is proximate to a quantum melting
point. Figure 1 schematically depicts the possible location of Ba3CoSb2O9 with
respect to a quantum critical point in a magnetic phase diagram.

Subsequent theoretical studies, including numerical tensor network
calculations [27] and $1/\mathcal{N}$ Schwinger boson (SB) calculations [28],
support the $XXZ$ model. Furthermore, the SB theory provides a novel
interpretation of the spectrum, suggesting that the low-lying energy magnetic
excitations consist of two spinon bound states, glued by emergent (Hubbard-
Stratonovich) gauge fields. These coexist with two high-energy dispersing
continua made of two quasi-free spinons. Namely, at high energies, the
confinement length of spinons is greater than unit space lattice. In the light
of these findings, it can be concluded that in the proximity of a quantum
melting point, the spin-1 collective excitation of the 120∘ Néel order seems
better described by a bound state of two spin-$\frac{1}{2}$ spinons —the
quasiparticles of the neighboring spin liquid— than the usual magnonic
excitations [29, 30].

Figure 1: Schematic phase diagram for quasi-2D frustrated antiferromagnets
[31, 32]. The abscissa axis ’g’ represents a generic measure of zero-point
quantum fluctuations, with ’gc’ connecting continuously the 120° Néel state to
a quantum spin liquid state at $T=$ 0\. Solid line denote the transition from
the classical renormalized to the quantum critical regime, marked by the Néel
temperature ($T_{N}$). The dashed line indicates the crossover around $T^{*}$
from the classical renormalized regime —with well-defined two-spinon bound
states across the Brillouin zone— to an intermediate regime where bound states
begin to break up. This region, termed the terminated Goldstone regime (TGR)
retains only bound states around the Goldstone modes. The blue and green
vertical bands depict the estimated positions of named compounds in the phase
diagram. The solid red line represents the location of the model used in our
calculations (see Section II) and the yellow region, emerging from gc, denote
a quantum critical regime.

More recently [33], the delafossite material KYbSe2, which realizes the
$J_{1}-J_{2}$ model on the triangular lattice, has been shown to be proximate
to a quantum critical point. In this study, the magnetic spectrum was observed
with INS at $T=$ 300 mK, slightly above the Néel temperature $T_{N}$ = 290 mK.
The spectrum displayed an extended continuum with a sharp lower edge that is
gapless at the momentum $K$, corresponding to the 120∘ Néel order expected
below $T_{N}$. The dynamical structure factor at the magnetic wave vector $K$
exhibits a scaling collapse in $~{}\hbar\omega/k_{B}T$ down to 0.3K,
indicative of a second-order quantum phase transition. The agreement with
theoretical predictions based on tensor network and $1/\mathcal{N}$ SB
calculations suggests that KYbSe2 is even closer to a quantum critical point
than Ba3CoSb2O9 (see Fig. 1). Based on the zero temperature SB prediction, it
is expected that the low-energy edge of the spectrum observed at 0.3K (greater
than $T_{N}$) will separate from the continuum as temperature passes through
the Néel temperature down, giving rise to collective excitations made of two
spinon bound states. However, it remains unclear whether this low-lying energy
band will be well-defined across the entire Brillouin zone. The consistency
found between theory and experiments in KYbSe2, along with entanglement
witnesses, strongly supports the strategy of investigating the magnetic
spectrum of compounds in proximity to a continuous quantum phase transition.
In this case, the indications are that KYbSe2 could be proximate to a Z2 spin
liquid state.

So far, the theoretical studies mentioned above have primarily been performed
at zero temperature, while only a few studies consider finite temperature [34,
35, 36]. To provide a more comprehensive description aligning with
experimental conditions, the next step is to incorporate the combined effects
of both interlayer exchange coupling and temperature. In compounds like
Ba3CoSb2O9 and KYbSe2, although the interlayer interaction is significantly
weaker than the intralayer coupling, it plays a crucial role in establishing a
finite Néel temperature.

In this work, we investigate the temperature effect on the magnetic excitation
spectra of spin 1/2 quasi-2D triangular AF. Unlike spin wave theory, the
Schwinger boson theory (SBT) relies on magnetic link degrees of freedom that
are believed to take better into account the mixing of transverse and
longitudinal fluctuations, as soon as Néel temperature is approached. Firstly,
we find that at the saddle point (SP) level the SBT recovers qualitatively
well the expected (static) behavior of the $120^{\circ}$ magnetically ordered
regime below Néel temperature. Secondly, at the $1/\mathcal{N}$ approximation
level, we have found an intriguing temperature dependence of the spectrum.
Specifically, we identify a low-temperature regime, below a crossover
temperature $T^{*}$, where the two-spinon bound state band is well-defined
across an ample region of the Brillouin zone. However, as temperature rises
above $T^{*}$, the gauge field fluctuations fail to form the two-spinon bound
states in a significant region of the Brillouin zone. This leads to a
differentiation in the magnetic excitation spectrum. Remarkably, for
temperatures around $T_{N}$, only two spinon bound states, centered around the
Goldstone modes, remain well-defined. This behavior is similar to what happens
in weakly 3D interacting antiferromagnetic chains [37]. In these systems,
clear signatures of quantum criticality are observed in the renormalized
classical regime (below but close to the Néel temperature). In analogy to the
superfluid 4He case, where phonon excitations terminate at some momentum and
energy by decaying into two roton excitations, we call this antiferromagnetic
region –between $T^{*}$ and $T_{N}$– terminated Goldstone regime (TGR) (see
Fig. 1), although the decay mechanisms are different [38]. These finite
temperature results strongly support the proposed behavior of frustrated
antiferromagnets near a quantum critical regime by Chubukov, Sachdev, and
Senthil [31], and provide a deeper insight into the fractionalization process
[39]

In Section II, we briefly developed the $1/\mathcal{N}$ Schwinger boson
formalism. Section III includes the interlayer coupling, enabling the study of
the Néel temperature. In Section IV, we analyze the effect of temperature on
the magnetic excitation spectrum. Section V concludes with final remarks.

##  II Schwinger boson formalism

In this section we extend the $1/\mathcal{N}$ SB theory to finite temperatures
for quasi two dimensional systems. In order to get a self contained
description we present the main steps of the calculation, emphasizing the
temperature dependence of the results. The details of the full calculation can
be found in Ref.[40, 41, 42, 30, 43].

We study the Heisenberg Hamiltonian in the layered triangular lattice

$\hat{H}=\sum_{<ij>}J_{ij}\hat{\textbf{S}}_{i}\cdot\hat{\textbf{S}}_{j},$ (1)

where $<ij>$ sums over all the nearest-neighbor pairs of spins of the layered
triangular lattice. The in-plane exchange interaction is $J$ while there is a
no frustrating interlayer exchange interaction $J_{\perp}$. For simplicity,
here we study the isotropic Heisenberg model with only nearest-neighbor
couplings. However, the theory can be straightforwardly extended to include
anisotropic [28] and further exchange interactions [33]. In fact, in the next
Section, we will present some Néel temperature estimations for real compounds
with more evolved Hamiltonian than the isotropic Heisenberg model.

Within the SB formalism the Heisenberg Hamiltonian (1) can be written exactly
as

$\hat{H}=\sum_{<ij>}J_{ij}\left(:\hat{B}^{\dagger}_{ij}\hat{B}_{ij}:-\hat{A}^{\dagger}_{ij}\hat{A}_{ij}\right),$
(2)

where
$\hat{A}_{ij}=\frac{1}{2}(\hat{b}_{i\uparrow}\hat{b}_{j\downarrow}-\hat{b}_{i\downarrow}\hat{b}_{j\uparrow})$
and
$\hat{B}_{ij}=\frac{1}{2}(\hat{b}_{i\uparrow}\hat{b}^{\dagger}_{j\downarrow}+\hat{b}_{i\downarrow}\hat{b}^{\dagger}_{j\uparrow})$
are the SU(2) invariant link operators of the theory. $A^{\dagger}_{ij}$
creates singlets between $ij$ sites and $B^{\dagger}_{ij}$ makes them
resonate. Furthermore, up to a constant,
$:\hat{B}^{\dagger}_{ij}\hat{B}_{ij}:\sim(\hat{\textbf{S}}_{i}+\hat{\textbf{S}}_{j})^{2}$
and
$\hat{A}^{\dagger}_{ij}\hat{A}_{ij}\sim(\hat{\textbf{S}}_{i}-\hat{\textbf{S}}_{j})^{2}$
which gives a measure of ferromagnetic and antiferromagnetic correlations,
respectively. Namely, the particular structure of Hamiltonian (2) is ideal to
treat Heisenberg models in the presence of magnetic frustration [44, 45, 46,
30].

The partition function for the model is expressed as the following path
integral over the coherent states of the SB’s

$\begin{split}\mathcal{Z}[j]=\int\mathcal{D}[\bar{b},b]\mathcal{D}[\lambda]&e^{-\int_{0}^{\beta}d\tau\left[\sum_{i\sigma}\bar{b}_{i\sigma}^{\tau}\partial_{\tau}{b}_{i\sigma}+\mathcal{H}(\bar{b},b)+\mathcal{J}_{b}+\mathcal{J}_{s}\right]}\\\
\times\ &e^{-\int_{0}^{\beta}d\tau\
i\sum_{i}\lambda_{i}^{\tau}\left(\sum_{\sigma}\bar{b}_{i\sigma}^{\tau}b_{i\sigma}-2S\right)},\end{split}$
(3)

where $\mathcal{H}(\bar{b},b)$ is the Hamiltonian evaluated at the site $i$
and time $\tau$ dependent complex eigenvalues of the coherent states,
$\mathcal{J}_{b}=\frac{1}{2}\sum_{i}h_{i}^{\mu}\textbf{b}_{i}^{\tau\dagger}\cdot\sigma^{\mu}\cdot\textbf{b}_{i}^{\tau}$
–with
$\textbf{b}^{\tau\dagger}_{i}=(\overline{b}^{\tau}_{i\uparrow},\overline{b}^{\tau}_{i\downarrow})$–
is the coupling between the spins and an infinitesimal symmetry breaking field
$\textbf{h}_{i}=h(\cos(\textbf{Q}\cdot\textbf{r}_{i}),\sin(\textbf{Q}\cdot\textbf{r}_{i}),0)$
that selects the $120^{\circ}$ Néel order in the x-y plane
[$\textbf{Q}=\left(\frac{2}{3}\pi,\frac{2}{\sqrt{3}}\pi,\pi\right)$] ,
$\mathcal{J}_{s}=\frac{1}{2}\sum_{i}j_{i}^{\tau\mu}\
\textbf{b}_{i}^{\tau\dagger}\cdot\sigma^{\mu}\cdot\textbf{b}_{i}^{\tau}$ is
the coupling between the spins and a external magnetic field $j_{i}^{\tau\mu}$
used to compute the dynamical spin susceptibility, and $\lambda_{i}^{\tau}$
are the Lagrange multipliers introduced to satisfy the local constraints
$b^{\dagger}_{i\uparrow}b_{j\uparrow}+b^{\dagger}_{i\downarrow}b_{j\downarrow}=2S$.

The standard procedure involves i) a Hubbard-Stratonovich (HS) transformation
to decouple the $A^{\dagger}_{ij}\\!A_{ij}$ and $B^{\dagger}_{ij}\\!B_{ij}$
terms (quartic in bosons) of $\mathcal{H}(\bar{b},b)$, ii) use a local
reference quantization axis for the spinors $b^{\tau}_{i\sigma}$ oriented in
the direction of the $120^{\circ}$ magnetic order Q, iii) integrate out the
quadratic bosonic part giving rise to

$\displaystyle\mathcal{Z}_{\rm bos}(\overline{W},W,\lambda,j)\\!\\!$
$\displaystyle=$
$\displaystyle\\!\\!\\!\int\\!D[\overline{b},b]e^{-\vec{b}^{\dagger}\cdot\mathcal{G}^{-1}(\overline{W},W,\lambda,j)\cdot\vec{b}}=$
$\displaystyle=$
$\displaystyle\det\left[\mathcal{G}(\overline{W},W,\lambda,j)\right],$

where $\vec{b}$ is a four component vector containing the variables
$b_{i\sigma}^{\tau}$ and $\overline{b}_{i\sigma}^{\tau}$, $\mathcal{G}^{-1}$
is the $4\times 4$ dynamical matrix (see next subsection), and
$W_{ij}^{\xi\tau}$ (with $\xi=A,\ B$) are the HS auxiliary fields, and finally
iv) exponentiate $\mathcal{Z}_{\rm bos}(\overline{W},W,\lambda,j)$, rendering
the total partition function as

$\mathcal{Z}[j]=\int\mathcal{D}[\overline{W},W]\mathcal{D}[\lambda]e^{-S_{\rm
eff}(\overline{W},W,\lambda,h,j)}$ (4)

with

${\textstyle S_{\rm
eff}(\overline{W},W,\lambda,h,j)=S_{0}(\overline{W},W,\lambda)+S_{\rm
bos}(\overline{W},W,\lambda,h,j)}$

where

$S_{0}(\overline{W},W,\lambda)=\int_{0}^{\beta}d\tau\big{(}\sum_{ij,\xi}J_{ij}\overline{W}_{ij}^{\xi\tau}W_{ij}^{\xi\tau}-i2S\sum_{i}\lambda_{i}^{\tau}\big{)}$
(5)

and

$\displaystyle S_{\rm bos}(\overline{W},W,\lambda,j)$ $\displaystyle\\!=\\!$
$\displaystyle-\frac{1}{2}\ln\ \mathcal{Z}_{\rm
bos}(\overline{W},W,\lambda,j)=$ (6) $\displaystyle\\!=\\!$
$\displaystyle\frac{1}{2}{\rm
Tr}\ln\left[\mathcal{G}^{-1}(\overline{W},W,\lambda,j)\right],$

with the trace, in the last line, taken over space, time, and boson flavor
indices. Equation (4) means that the interacting spin problem has been mapped
exactly to a problem of free spin-$\frac{1}{2}$ bosons interacting with
fluctuating (space and time) auxiliary fields $W_{ij}^{A\tau}$,
$W_{ij}^{B\tau}$, and $\lambda_{i}^{\tau}$ which are the gauge fields of the
theory [30].

Next, a saddle point expansion is performed and, keeping up to the Gaussian
fluctuations of the auxiliary fields, the partition function yields

$\mathcal{Z}^{(2)}\\![j]\\!=\\!e^{-S^{\rm sp}_{\rm eff}{(\overline{W}_{\rm
sp},{W}_{\rm sp},\lambda_{\rm
sp},j)}}\\!\times\\!\int\\!\\!D[\bar{\phi},\phi]\\!\
e^{-\Delta\vec{\phi}^{\dagger}\cdot S^{(2)}\cdot\Delta\vec{\phi}},$ (7)

where the first exponential is the partition function within the SP
approximation, and the integral corresponds to the contribution of the
Gaussian fluctuations of the auxiliary fields
$\Delta\vec{\phi}^{\dagger}=\vec{\phi}^{\dagger}-\vec{\phi}^{{\rm sp}\dagger}$
with $\vec{\phi}^{\dagger}=\left(\overline{W}_{i,\delta}^{\xi\tau},\
W_{i,\delta}^{\xi\tau},\ \lambda_{i}^{\tau}\right)$. $S^{(2)}$ is the
fluctuation matrix composed by the second derivatives of the effective action
$S_{\rm eff}$ with respect to the auxiliary fields [30]. It is worth to stress
that once the SB theory is extended to a large $\mathcal{N}$ number of
flavors, the contribution of Gaussian fluctuations to the dynamical spin
susceptibility is of order $1/\mathcal{N}$, so that, we call the present
theory $1/\mathcal{N}$ SB theory although we are working within the physical
case $\mathcal{N}=2$ [40, 41, 43].

###  II.1 Saddle point equations

The SP approximation requires that

$\frac{\partial S_{eff}}{\partial\phi_{\alpha}}\bigg{|}_{sp}=\frac{\partial
S_{0}}{\partial\phi_{\alpha}}\bigg{|}_{sp}+\frac{1}{2}\text{Tr}\left[\mathcal{G}_{sp}\frac{\partial\mathcal{G}^{-1}}{\partial\phi_{\alpha}}\bigg{|}_{sp}\right]=0.$
(8)

where $\phi_{\alpha}$ represents a given HS field (including $\lambda$) at
some position and time. Taking the ansatz compatible with the magnetic order
in the $xy$ plane $\langle
A_{\bm{\delta}}\rangle=-\langle\overline{A}_{\bm{\delta}}\rangle=iA_{\bm{\delta}}$
, $\langle
B_{\bm{\delta}}\rangle=\langle\overline{B}_{\bm{\delta}}\rangle=B_{\bm{\delta}}$
($\bm{\delta}$ represents the vector that connects sites $i$ and $j$) and
$\lambda_{sp}=i\lambda$, the dynamical matrix in the frequency and momentum
representation, using the spinor basis
$\textbf{b}^{\omega\dagger}_{\textbf{k}}=(\bar{b}^{\omega}_{\textbf{k}\uparrow},{b}^{-\omega}_{-\textbf{k}\downarrow},\bar{b}^{\omega}_{\textbf{k}\downarrow},{b}^{-\omega}_{-\textbf{k}\uparrow})$,
is written as

$\scriptsize\mathcal{G}^{-1}_{sp}(\textbf{k},i\omega)=\begin{pmatrix}\left(i\omega+\lambda_{sp}+\gamma^{B}_{\textbf{k}+\frac{\textbf{Q}}{2}}\right)e^{-i\omega
0^{+}}&-\gamma^{A}_{\textbf{k}+\frac{\textbf{Q}}{2}}&\frac{h}{2}&0\\\
-\gamma^{A}_{\textbf{k}+\frac{\textbf{Q}}{2}}&\left(-i\omega+\lambda_{sp}+\gamma^{B}_{\textbf{k}+\frac{\textbf{Q}}{2}}\right)e^{i\omega
0^{+}}&0&\frac{h}{2}\\\
\frac{h}{2}&0&\left(i\omega+\lambda_{sp}+\gamma^{B}_{-\textbf{k}+\frac{\textbf{Q}}{2}}\right)e^{-i\omega
0^{+}}&-\gamma^{A}_{-\textbf{k}+\frac{\textbf{Q}}{2}}\\\
0&\frac{h}{2}&-\gamma^{A}_{-\textbf{k}+\frac{\textbf{Q}}{2}}&\left(-i\omega+\lambda_{sp}+\gamma^{B}_{-\textbf{k}+\frac{\textbf{Q}}{2}}\right)e^{i\omega
0^{+}}\end{pmatrix}$ (9)

with functions
$\gamma_{\bm{k}}^{A}\\!\\!\\!=\\!\\!\\!\sum_{\bm{\delta}>0}J_{\bm{\delta}}A_{\bm{\delta}}\sin(\textbf{k}\cdot\bm{\delta})$
and
$\gamma_{\bm{k}}^{B}=\sum_{\bm{\delta}>0}J_{\bm{\delta}}B_{\bm{\delta}}\cos(\textbf{k}\cdot\bm{\delta})$.
Convergence factors $e^{\pm i\omega 0^{+}}$ are necessary to correctly take
into account the temporal ordering within the functional integral. The inverse
of the dynamical matrix at the SP is the spinon propagator $\mathcal{G}_{sp}$
that can be decomposed in term of simple fractions as

$\mathcal{G}_{sp}(\textbf{k},i\omega)=\sum_{\sigma\sigma^{\prime}}\frac{g^{\sigma\sigma^{\prime}}\scriptstyle{(\textbf{k})}}{i\omega+\sigma\varepsilon_{\textbf{k}}^{\sigma^{\prime}}(h)}$
(10)

whose poles are the dispersion relations of the free spinons
$\varepsilon_{\textbf{k}}^{\sigma}(h)$, with $\sigma,\sigma^{\prime}=\pm$,
that can be obtained together with the matrices $g^{\sigma\sigma^{\prime}}$
after performing a paraunitary diagonalization of the dynamical matrix (See
Appendix A). Replacing the expressions of $\mathcal{G}^{-1}_{sp}$ and
$\mathcal{G}_{sp}$ in the SP equation (8), and carrying out the Matsubara
sums, we arrive at the self-consistent set of equations

$\displaystyle A_{\bm{\delta}}$
$\displaystyle=\frac{1}{2N_{s}}\sum_{\textbf{k}}\tilde{A}_{\textbf{k}}\sin(\textbf{k}+\frac{\textbf{Q}}{2})\cdot\bm{\delta}$
(11) $\displaystyle B_{\bm{\delta}}$
$\displaystyle=\frac{1}{2N_{s}}\sum_{\textbf{k}}\tilde{B}_{\textbf{k}}\cos(\textbf{k}+\frac{\textbf{Q}}{2})\cdot\bm{\delta}$
(12) $\displaystyle S$
$\displaystyle=\frac{1}{2N_{s}}\sum_{\textbf{k}}\tilde{B}_{\textbf{k}}$ (13)

where $N_{s}$ is the total number of sites and

$\begin{split}\tilde{A}_{\textbf{k}}&=\sum_{\sigma\sigma^{\prime}}\left(g^{\sigma\sigma^{\prime}}_{21}{\scriptstyle{(\textbf{k})}}+g^{\sigma\sigma^{\prime}}_{43}\scriptstyle{(-\textbf{k})}\right)\left[1+n(\sigma\varepsilon_{\textbf{k}}^{\sigma^{\prime}}(h))\right]\\\
\tilde{B}_{\textbf{k}}&=\sum_{\sigma\sigma^{\prime}}g^{\sigma\sigma^{\prime}}_{11}{\scriptstyle{(\textbf{k})}}\
n(\sigma\varepsilon_{\textbf{k}}^{\sigma^{\prime}}(h))+g^{\sigma\sigma^{\prime}}_{22}{\scriptstyle(\textbf{k})}\left[1+n(\sigma\varepsilon_{\textbf{k}}^{\sigma^{\prime}}(h))\right].\end{split}$
(14)

The explicit temperature dependence is present in the Bose occupation numbers
$n$.

###  II.2 Dynamical spin susceptibility

The dynamical spin susceptibility is obtained as [41, 30]

$\chi_{\mu\nu}(\textbf{q},i\omega)=\lim_{h\rightarrow
0}\lim_{N_{s}\rightarrow\infty}\frac{\partial^{2}\ln\mathcal{Z}[j]}{\partial
j^{\mu}_{\textbf{q},i\omega}\partial
j^{\nu}_{-\textbf{q},-i\omega}}\bigg{|}_{j=0},$ (15)

where the order of the limits is a key point of the calculation since we are
interested in the excitation spectrum of magnetically ordered states proximate
to a quantum melting point. In particular, to compute the dynamical spin
susceptibility up to Gaussian order –above the SP solution–
$\mathcal{Z}^{(2)}\\![j]$ of Eq. (7) is plugged in Eq. (15), giving rise to
three terms of order $1/\mathcal{N}$. In Refs. [30, 43] we have shown that by
only keeping the bubble like diagram the dynamical susceptibility can be
written as

$\chi_{\mu\mu}(\textbf{q},i\omega)=\chi^{sp}_{\mu\mu}(\textbf{q},i\omega)+\chi^{fl}_{\mu\mu}(\textbf{q},i\omega),$
(16)

being

$\chi^{sp}_{\mu\mu}(\textbf{q},i\omega)=\frac{1}{2}\text{Tr}\left[\mathcal{G}_{sp}\frac{\partial\mathcal{G}^{-1}}{\partial
j^{u}_{\textbf{q},i\omega}}\bigg{|}_{sp}\mathcal{G}_{sp}\frac{\partial\mathcal{G}^{-1}}{\partial
j^{u}_{-\textbf{q},-i\omega}}\bigg{|}_{sp}\right]$ (17)

and

$\chi^{fl}_{\mu\mu}(\textbf{q},i\omega)=\sum_{\alpha_{1}\alpha_{2}}\Lambda^{\mu}_{\phi_{\alpha_{1}}}(\textbf{q},i\omega)\mathcal{D}_{\alpha_{2}\alpha_{1}}(\textbf{q},i\omega)\Lambda^{\mu}_{\phi_{\alpha_{2}}}(-\textbf{q},-i\omega)$
(18)

where

$\Lambda^{\mu}_{\phi_{\alpha_{i}}}(\textbf{q},i\omega)=\frac{1}{2}\text{Tr}\left[\mathcal{G}_{sp}\frac{\partial\mathcal{G}^{-1}}{\partial\phi_{\alpha_{i}}}\bigg{|}_{sp}\mathcal{G}_{sp}\frac{\partial\mathcal{G}^{-1}}{\partial
j^{u}_{\textbf{q},i\omega}}\bigg{|}_{sp}\right],$ (19)

and $\mathcal{D}_{\alpha_{1}\alpha_{2}}$ is the inverse of the fluctuation
matrix
$S^{(2)}_{\alpha_{1}\alpha_{2}}=\frac{1}{2}\frac{\partial^{2}S_{eff}}{\partial\phi_{\alpha_{1}}\partial\phi_{\alpha_{2}}}\big{|}_{sp}=\frac{1}{2}\left\\{\frac{\partial^{2}S_{0}}{\partial\phi_{\alpha_{1}}\partial\phi_{\alpha_{2}}}-\Pi_{\alpha_{1}\alpha_{2}}\right\\}$
with

$\Pi_{\alpha_{1}\alpha_{2}}=\frac{1}{2}\text{Tr}\left[\mathcal{G}_{sp}\frac{\partial\mathcal{G}^{-1}}{\partial\phi_{\alpha_{1}}}\bigg{|}_{sp}\mathcal{G}_{sp}\frac{\partial\mathcal{G}^{-1}}{\partial\phi_{\alpha_{2}}}\bigg{|}_{sp}\right].$
(20)

After carrying out the Matsubara’s frequency sums in (17), (19) and (20) for
the out of plane $zz$ entry of the susceptibility, we arrive at

$\begin{split}\chi_{zz}^{sp}(\textbf{q},i\omega)=\frac{1}{8N_{s}}\sum_{\textbf{k}}\sum_{\sigma\sigma^{\prime}\sigma^{\prime\prime}}\Bigg{\\{}&\frac{Tr\left[g^{\sigma\sigma^{\prime}}{\scriptstyle{(\textbf{k}+\textbf{q})}}u^{z}g^{\sigma\sigma^{\prime\prime}}{\scriptstyle{(\textbf{k})}}u^{z}\right]}{\varepsilon^{\sigma^{\prime}}_{\textbf{k}+\textbf{q}}-\varepsilon^{\sigma^{\prime\prime}}_{\textbf{k}}+\sigma
i\omega}\left(n(\varepsilon^{\sigma^{\prime\prime}}_{\textbf{k}})-n(\varepsilon^{\sigma^{\prime}}_{\textbf{k}+\textbf{q}})\right)-\\\
&\frac{Tr\left[g^{\sigma\sigma^{\prime}}{\scriptstyle{(\textbf{k}+\textbf{q})}}u^{z}g^{\bar{\sigma}\sigma^{\prime\prime}}{\scriptstyle{(\textbf{k})}}u^{z}\right]}{\varepsilon^{\sigma^{\prime}}_{\textbf{k}+\textbf{q}}+\varepsilon^{\sigma^{\prime\prime}}_{\textbf{k}}+\sigma
i\omega}\left(n(\varepsilon^{\sigma^{\prime\prime}}_{\textbf{k}})+n(\varepsilon^{\sigma^{\prime}}_{\textbf{k}+\textbf{q}})+1\right)\Bigg{\\}}\end{split}$
(21)
$\begin{split}\Lambda^{z}_{\phi_{\alpha_{1}}}(\textbf{q},i\omega)=\frac{1}{4N_{s}}\sum_{\textbf{k}}\sum_{\sigma\sigma^{\prime}\sigma^{\prime\prime}}\Bigg{\\{}&\frac{Tr\left[g^{\sigma\sigma^{\prime}}{\scriptstyle{(\textbf{k}+\textbf{q})}}u^{z}g^{\sigma\sigma^{\prime\prime}}{\scriptstyle{(\textbf{k})}}v^{\textbf{k},\textbf{k}+\textbf{q}}_{\phi_{\alpha_{1}}}\right]}{\varepsilon^{\sigma^{\prime}}_{\textbf{k}+\textbf{q}}-\varepsilon^{\sigma^{\prime\prime}}_{\textbf{k}}+\sigma
i\omega}\left(n(\varepsilon^{\sigma^{\prime\prime}}_{\textbf{k}})-n(\varepsilon^{\sigma^{\prime}}_{\textbf{k}+\textbf{q}})\right)-\\\
&\frac{Tr\left[g^{\sigma\sigma^{\prime}}{\scriptstyle{(\textbf{k}+\textbf{q})}}u^{z}g^{\bar{\sigma}\sigma^{\prime\prime}}{\scriptstyle{(\textbf{k})}}v^{\textbf{k},\textbf{k}+\textbf{q}}_{\phi_{\alpha_{1}}}\right]}{\varepsilon^{\sigma^{\prime}}_{\textbf{k}+\textbf{q}}+\varepsilon^{\sigma^{\prime\prime}}_{\textbf{k}}+\sigma
i\omega}\left(n(\varepsilon^{\sigma^{\prime\prime}}_{\textbf{k}})+n(\varepsilon^{\sigma^{\prime}}_{\textbf{k}+\textbf{q}})+1\right)\Bigg{\\}}\end{split}$
(22)
$\begin{split}\Pi_{\alpha_{1}\alpha_{2}}(\textbf{q},i\omega)=\frac{1}{2N_{s}}\sum_{\textbf{k}}\sum_{\sigma\sigma^{\prime}\sigma^{\prime\prime}}\Bigg{\\{}&\frac{Tr\left[g^{\sigma\sigma^{\prime}}{\scriptstyle{(\textbf{k}+\textbf{q})}}v^{\textbf{k}+\textbf{q},\textbf{k}}_{\phi_{\alpha_{1}}}g^{\sigma\sigma^{\prime\prime}}{\scriptstyle{(\textbf{k})}}v^{\textbf{k},\textbf{k}+\textbf{q}}_{\phi_{\alpha_{2}}}\right]}{\varepsilon^{\sigma^{\prime}}_{\textbf{k}+\textbf{q}}-\varepsilon^{\sigma^{\prime\prime}}_{\textbf{k}}+\sigma
i\omega}\left(n(\varepsilon^{\sigma^{\prime\prime}}_{\textbf{k}})-n(\varepsilon^{\sigma^{\prime}}_{\textbf{k}+\textbf{q}})\right)-\\\
&\frac{Tr\left[g^{\sigma\sigma^{\prime}}{\scriptstyle{(\textbf{k}+\textbf{q})}}v^{\textbf{k}+\textbf{q},\textbf{k}}_{\phi_{\alpha_{1}}}g^{\bar{\sigma}\sigma^{\prime\prime}}{\scriptstyle{(\textbf{k})}}v^{\textbf{k},\textbf{k}+\textbf{q}}_{\phi_{\alpha_{2}}}\right]}{\varepsilon^{\sigma^{\prime}}_{\textbf{k}+\textbf{q}}+\varepsilon^{\sigma^{\prime\prime}}_{\textbf{k}}+\sigma
i\omega}\left(n(\varepsilon^{\sigma^{\prime\prime}}_{\textbf{k}})+n(\varepsilon^{\sigma^{\prime}}_{\textbf{k}+\textbf{q}})+1\right)\Bigg{\\}}\end{split}$
(23)

where $u_{z}=diag(1,-1,-1,1)$ come from the external vertices
$\frac{\partial\mathcal{G}^{-1}}{\partial j^{u}_{\textbf{q},i\omega}}$ and
$v^{\textbf{k},\textbf{k}^{\prime}}_{\phi_{\alpha}}$ come from the internal
vertices $\frac{\partial\mathcal{G}^{-1}}{\partial\phi_{\alpha}}$. Similar
expressions can be obtained for the in-plane ($xx$ and $yy$) entries of the
susceptibility.

Finally, by means of the fluctuation-dissipation theorem, we can obtain the
dynamical spin structure factor by performing the analytic continuation
$i\omega\rightarrow\omega+i\eta$ in the dynamical spin susceptibility

$S^{\mu\mu}(\textbf{q},\omega)=\left(\frac{-1}{\pi}\right)\frac{1}{1-e^{-\beta\omega}}\text{Im}\big{[}\chi^{\mu\mu}(\textbf{q},\omega+i\eta)\big{]},$
(24)

where $\eta$ is the artificial broadening.

##  III Magnetically ordered regime

In this section we focus on the ability of the SBT to describe the finite
temperature behavior as a function of $J_{\perp}$. In particular, in the
ordered regime, we find that the SP solution recovers qualitatively well the
expected (static) behavior of quasi-2D triangular AFs.

In accordance with the Mermin-Wagner theorem [40, 41], the onset of long-range
magnetic order in a triangular Heisenberg antiferromagnet –with continuous
symmetry and short range exchange interactions– can only occur at zero
temperature. Correspondingly, by mapping the 2D Heisenberg antiferromagnet to
the effective non linear sigma model (NLSM), it is shown that there is a
crossover from a high temperature paramagnetic regime to a low temperature
classical renormalized regime where the finite temperature spin correlation
functions increases exponentially as temperature goes down to zero [47]. The
SBT reproduces quite well this low temperature behavior of the spin
correlations at the saddle point level [48, 31]. In particular, it recovers
the low temperature $120^{\circ}$ Néel correlation functions predicted by
NLSM, while at zero temperature the condensation of the SBs signals the
occurrence of a finite magnetization. At higher temperatures, however, the SP
approximation overestimates the entropy –due to violation of the local
constraint– so neighboring spins become perfectly uncorrelated at certain
temperature [49]. This unphysical result has been studied in the context of
large $\mathcal{N}$ theory along with the possibility of being corrected for
finite $\mathcal{N}$ [50].

On the other hand, any finite interlayer exchange $J_{\perp}$ gives rise to
long range magnetic order below a finite Néel temperature. The very definition
of a quasi-2D antiferromagnet implies that $J_{\perp}$ is much smaller than
the in-plane exchange interaction $J$ –two or more orders of magnitude
smaller–, making it difficult to model the compounds in order to extract
$J_{\perp}$ values from the measured Néel temperatures.

Formally, the magnetic phase transition is derived by studying the behavior of
the SB self consistent Eqs. (11-13) in the thermodynamic limit, namely, first
taking $N_{s}\rightarrow\infty$ and then $h\rightarrow 0$ [41]. In particular,
there is a regime where $\varepsilon^{+}_{{\bf 0}}\rightarrow 0$ and
$\varepsilon^{-}_{{\bf 0}}\neq 0$ which corresponds to the condensation of the
SB $\sigma=+$ at ${\bf q}=0$. This condensation is directly related to the
presence of a local magnetization that lies in the $xy$ plane and points in
the $x$ direction of a local reference frame that rotates 120∘ from site to
site. In practice, this calculations is carried out on finite size systems of
$N_{s}=N_{\perp}\times N_{\parallel}$ sites, where $N_{\perp}$ is the number
of layers taken and $N_{\parallel}$ is the number of sites in each triangular
lattice layer. Then, a symmetry breaking field $h=1/N_{\parallel}$ is applied.
Since $\varepsilon^{+}_{{\bf 0}}\sim 1/N_{\parallel}$ and
$\varepsilon^{-}_{{\bf 0}}\sim(1/N_{\parallel})^{\frac{1}{2}}$, the condensate
is identified with the ${\bf q}=0$ and $+$ term of the constraint equation
(13). Fig. 2 shows the dependence of this condensate with temperature for a
quasi-2D triangular system composed by 8 layers –with periodic boundary
conditions– and $J_{\perp}/J=0.01$. In this work, we have taken
$N_{\parallel}=144\times 144$. This figure also shows the spinon relation
dispersion values $\varepsilon^{+}_{\bf 0}$ and $\varepsilon^{-}_{{\bf 0}}$
where it can be seen that $\varepsilon^{+}_{{\bf 0}}\sim 0$ (gapless) and
$\varepsilon^{-}_{{\bf 0}}\neq 0$ (gapped) for $T\leq T_{N}$ (consistent with
$\sigma=+$, $\bm{q}=\bm{0}$ condensate); while for $T>T_{N}$, it is opened a
gap in both flavors $\varepsilon^{+}_{{\bf 0}}=\varepsilon^{-}_{{\bf 0}}\neq
0$ what signals no condensation, that is, the absence of magnetic order. In
the next section we will see that this behavior of the condensate has a
crucial effect on the magnetic spectrum as one approaches to the Néel
temperature.

Figure 2: $A_{\perp}$ (blue circles), condensate (purple triangles),
$\sigma\\!=\\!+$ (red diamonds) and $\sigma\\!=\\!-$ (black squares) spinon
dispersion relation gaps as a function of the temperature for $N_{s}=144\times
144\times 8$ sites and $J_{\perp}=0.01J$. Inset: Temperature dependence of the
in-plane saddle point parameters $A_{\parallel}$ and $B_{\parallel}$.

One way to measure the spin correlation between nearest-neighbor spins along
the interlayer direction is by computing the saddle-point value of the
Hubbard-Stratonovich field $A_{\perp}$, through the relation

$\langle\bm{S}_{i}\cdot\bm{S}_{i+c}\rangle=B_{\perp}^{2}-A_{\perp}^{2}.$ (25)

As the interlayer exchange interactions are not frustrated, $B_{\perp}$
vanishes, signalling the antiparallel spin alignment along the perpendicular
direction. The behavior of $A_{\perp}$ versus temperature is shown in Fig. 2.
There we observe that $A_{\perp}$ behaves similar to the condensate,
particularly vanishing at the Néel temperature, so that the fact that
$A_{\perp}\neq 0$ is an evidence of the 3D character of the magnetic order.
Namely, if $A_{\perp}\neq 0$, the system exhibits long-range magnetic order,
whereas above $T_{N}$ the triangular layers are magnetically decoupled
($A_{\perp}=0$) and the local magnetization goes to zero, in accordance with
Mermin-Wagner physics. In other words, the transition between the magnetically
ordered and paramagnetic phase occurs due to the effective
bidimensionalization of the magnetic state as temperature increases. So, the
vanishing of $A_{\perp}$ is another good signal of the Néel temperature.
Notice that as the temperature increases, the spin correlations in all
directions are weakened. One would expect, naively, that $A_{\perp}$ is of the
order of $J_{\perp},$ however any finite value of the interlayer coupling is
enough to generate strong correlations between triangular layers
($A_{\perp}\simeq 0.3$ at $T=0$).

The inset of Fig. 2 shows the dependence of the in-plane parameters
$A_{\parallel}$ and $B_{\parallel}$ with temperature. Even if the SP solution
predicts the above mentioned unphysical transition to a perfectly uncorrelated
paramagnet around $T\sim 0.43J$, when $T\lesssim 0.2J$ the SP solution
predicts reasonable results for the magnetically ordered regime of quasi-2D
triangular AFs.

Figure 3: Néel temperature as a function of the interlayer coupling
$J_{\perp}$ for $N_{\parallel}=144$ and different number of layers $N_{\perp}$
(with periodic boundary conditions). Inset: The solid red line corresponds to
the non-linear fit (26) for the 8 layers SB result.

Fig. 3 shows the Néel temperature as a function of $J_{\perp}$ and the number
of triangular layers with periodic boundary conditions. For very low
$J_{\perp}$ it is observed that $T_{N}$ remains practically independent of the
number of layers for $N_{\perp}\geq 4$, simplifying notoriously the
calculation of $T_{N}$ for quasi-2D systems. It is worth to stress that linear
spin wave theory applied to the quasi-2D triangular AF predicts larger values
of the Néel temperatures than the SBT ones [51]. For instance, for
$J_{\perp}/J=0.25$, $T^{LSWT}_{N}/J\sim 0.8$ while $T^{SBT}_{N}/J\sim 0.32$.
This clearly indicates that LSWT strongly overestimates the finite temperature
magnetization with respect to SBT. These values are consistent with the
expected behavior near $T_{N}$, where the transverse and longitudinal
fluctuations become of the same order. This feature is surely better captured
by the link fields $A_{ij}$ and $B_{ij}$ of the SBT than the local description
of LSWT. In addition, the SB results are compared with a non linear fit
corresponding to the empirical formula

$T_{N}(J_{\perp})=\frac{4\pi\rho_{s}}{b-\ln(J_{\perp})}$ (26)

based on the random phase approximation combined with quantum Monte Carlo
simulations performed in unfrustrated Heisenberg AFs [52]. Here, the best fit
is given by $\rho_{s}=0.12675$, $b=3.33257$. Notably, the SB results follows
quite well the empirical formula down to $J_{\perp}\sim 10^{-3}$ (see inset of
Fig. 3). Taking into account the complication of perfoming QMC calculations at
very low temperatures on triangular lattices –due to sign problem– the fact
that the best fit gives a value for the spin stiffness that is close to the SP
result already computed at zero temperature, $\rho_{s}=0.11$ [53], lends
support to our procedure to estimate $J_{\perp}$ values from $T_{N}$.

Figure 4: Total dynamical structure factor (DSF) for the quasi-2D triangular
antiferromagnet ($J_{\perp}=0.01J$) obtained with the Schwinger boson theory
at different temperatures. Left column: saddle point contribution to the DSF.
Right column: saddle point plus Gaussian fluctuations contribution to the DSF.
At the Néel temperature $T=0.2J$, Gaussian fluctuations vanish and the
remaining signal comes from the saddle point. Notice that there is a small gap
at $\Gamma$ as the spectra are computed for $k_{z}=\pi$. The open circles at
the top right panel aid to identify the two-spinon bound state band,
corresponding to the $S^{zz}$ component. The artificial broadening is taken
$\eta=0.02$.

The above results for the Heisenberg model can be easily extended to the $XXZ$
and $J_{1}-J_{2}$ models. In both cases we have observed that TN behaves in
the expected manner, that is, for a fixed value of $J_{\perp}$, TN increases
as anisotropy $\Delta$ increases; whereas TN decreases as frustration
$J_{2}/J_{1}$ increases. In Ba3CoSb2O9, described by the $XXZ$ model with
estimated values $J^{z}/J=0.937$ and $J_{\perp}/J=0.061$, the Néel temperature
gives $T^{\rm XXZ}_{N}/J=0.274$, which should be compared with the
experimental value $T^{\rm exp}_{N}/J^{\rm exp}\sim 0.2$ [25, 26, 28]. In
KYbSe2, described by the $J_{1}-J_{2}$ model with estimated value
$J_{2}/J_{1}=0.05$ and the Néel temperature $T^{\rm exp}_{N}/J_{1}^{\rm
exp}\sim 0.045$ [33] we obtain a magnitude for $J_{\perp}$ that is smaller
than $0.01\%$ of $J$. This result is in line with the experimental certainty
that KYbSe2 is one of the best examples of a two dimensional triangular AF.

Finally, it is important to remark that although we do not expect to obtain a
precise prediction for the values of $J_{\perp}$ or $T_{N}$ from a saddle
point (mean field) approximation, the fact that the SP solution gives the
expected qualitative description of finite temperature properties is crucial
to compute reliable spectra (to order $1/\mathcal{N}$) in order to get a
consistent thermal picture of the magnetic excitations.

##  IV Spectra

###  IV.1 Temperature evolution

In this subsection we focus on the temperature evolution of the magnetic
spectra. Before, we recall our previous results regarding the zero temperature
predictions of SBT for the excitation spectrum of the 2D triangular Heisenberg
model [46, 30, 43]. At the saddle point level [46] the excitation spectrum
consist of a two free spinon continuum (see top panel of first column of Fig
4). The lower edge of the continuum corresponds to
$\min(\varepsilon^{\sigma^{\prime}}_{\bm{k}+\bm{q}}+\varepsilon^{\sigma}_{\bm{k}})$
(see denominator in Eq. 21) and it is obtained when one spinon excitation is
created in the normal phase and another in the condensate. Furthermore, due to
the local constraint violation, the magnetic spectrum exhibit spurious modes
(see the gapped band at momentum $K$) arising from unphysical density
fluctuations [54]. The Gaussian correction changes drastically the SP spectrum
[30, 28]: i) it cancels out both, the lower edge of the SP two spinon
continuum and the spurious band and ii) it introduces new collective modes
(the poles of the RPA propagator $\mathcal{D}$) consisting of low energy two-
spinon bound states along with a quasi-free two spinon continuum at higher
energy (see top panel of second column of Fig 5). These two-spinon bound
states are the single-magnon excitations since, in the large $S$ limit, we
have numerically found that the $1/\mathcal{N}$ SBT spectrum coincides with
the LSWT prediction [43]. Namely, in the large $S$ limit, the two-spinon
continuum vanishes and the two-spinon bound state band coincides exactly with
both, the magnon band and the spectral weight predicted by LSWT. Then, for the
physical case of spin 1/2, it is fair to identify the two-spinon bound state
with the single-magnon excitations.

Fig. 4 displays the temperature evolution of the magnetic spectra along the
Brillouin zone path shown in the $T=0$ panel. On the left column, the saddle-
point predictions are shown, and it can be seen that, there is almost no
change in the structure and intensity of the spectra, except very close to the
Néel temperature (bottom left panel). In contrast, the
$1/\mathcal{N}$-corrected results, shown in the right column, exhibits a non-
trivial behavior with temperature. As temperature increases, a continuum
signal starts filling the gap between the single-magnon dispersion and the
$T=0$ continuum, due to the thermal creation and absorption of spinons (see
denominators
$\varepsilon^{\sigma^{\prime}}_{\bm{k}+\bm{q}}-\varepsilon^{\sigma}_{\bm{k}}$
in Eq. 23). These thermal processes take place at the same energy of the
magnon band, destroying the coherence of the two-spinon bound states, which
decay to two quasi-free spinon states. Consequently, the single-magnon
excitations acquire a finite lifetime, signalled by the broadening of the
magnon quasiparticle (see below). In a semiclassical picture, at $T=0$, due to
the non-collinear magnetic order, the transverse and longitudinal magnetic
fluctuations are tightly coupled. Then, the effect of finite temperature is to
incoherently decouple both kind fluctuations and the $1/\mathcal{N}$
corrections correctly capture this behavior.

We identify a crossover temperature $T^{*}\approx 0.75T_{N}$, above which the
magnon excitations loose their coherence across a significant region of the
first Brillouin zone, due to the mentioned decay in the emergent continuum.
Only in the vicinity of the Goldstone modes, magnons survive as well-defined
quasiparticles. Ultimately, above the Néel temperature, the spectra becomes
fully diffusive, and the $1/\mathcal{N}$ corrections to the magnetic dynamical
susceptibility vanish. Notice that saddle-point and $1/\mathcal{N}$-corrected
results are practically the same (see $T=0.20$ panel in Fig. 4). This is
equivalent to saying that, when temperature is high enough, the $SU(2)$
symmetry is restored.

Figure 5: Thermal evolution of the transverse dynamical structure factor
$S^{zz}(\textbf{q},\omega)$ for the quasi 2D triangular antiferromagnet
($J_{\perp}=0.01J$) obtained with the Schwinger boson theory. Inside the
dashed square it is shown how the emergent thermal continuum mixes with the
magnon band at $T=0.15J$, producing its decay.

Figure 5 shows the component $S^{zz}(\bm{q},\omega)$ of the dynamical spin
structure factor along the same path of the Brillouin zone as Fig. 4. As the
long-range magnetic order lies in the $x-y$ plane, $S^{zz}(\bm{q},\omega)$
takes only into account transverse magnetic fluctuations, whereas the total
spectra contains also the longitudinal fluctuations, that for non-collinear
orders is tightly coupled with the in-plane transverse ones. It can be noticed
that the relative spectral weight of the continuum signal with respect to the
single-magnon signal is smaller for the $zz$ component than for the total one
(see right column in Fig. 4). In other words, the anomalous extended continua,
characteristic of the triangular quantum antiferromagnet, has its origin
mainly in the in-plane spin fluctuations. The white dashed square is a visual
aid to help following the thermal evolution of a sector of the single-magnon
dispersion, between $K$ and $M$ points. For low temperatures, the magnon
excitations in this sector are well defined as there is almost no change with
respect to $T=0$. However, around $T^{*}\simeq 0.15$, the emergent thermal
spinons produce the decay of the two-spinon bound states (magnons). Beyond
$T^{*}$, the whole spectra become diffusive except in the neighborhood of the
Goldstone modes.

Figure 6: Thermal evolution of the transverse dynamical structure factor
$S^{zz}(\textbf{q},\omega)$ for (a)
$\textbf{q}_{1}\\!=\\!(\frac{3\pi}{4},\frac{7\pi}{4\sqrt{3}},\pi)$ and (b)
$\textbf{q}_{2}\\!=\\!(\frac{5\pi}{6},\frac{3\pi}{2\sqrt{3}},\pi)$. The two
momenta are shown in a plot of the BZ path at the top right corner of the (b)
panel. The magnon peak signal holds to higher temperatures for the closest
momentum to the Goldstone mode, $\textbf{q}_{1}$. At the Néel temperature
$T\\!=\\!0.2J$, the fluctuations vanish completely, and the excitations are
only due to free spinons.

###  IV.2 Relative spectral weight of the two-spinon bound states

To support our picture about the thermal evolution of the magnons and to
identify the crossover temperature, we analyze the spectral weight
contribution of the two-spinon bound states and fix a criterion to quantify
$T^{*}$. Fig. 6 displays $S^{zz}(\bm{q},\omega)$ for two different momenta
$\bm{q}_{1}$ and $\bm{q}_{2}$. It can be seen that for $\bm{q}_{1}$, close to
Goldstone point $K$, although the magnon loses spectral weight as the
temperature increases, it is always well defined, up to the Néel temperature
($T_{N}=0.2$). On the other hand, for $\bm{q}_{2}$ between the $K$ and $M$
points, the quasiparticle coherence is lost for $T\gtrsim 0.15$, being its
spectral weight transferred to the continuum background. This behaviour is
observed in an ample region of Brillouin zone. In order to see this, in Fig. 7
we compare the magnon spectral weight with the total signal of the dynamical
structure factor along the $K\to M$ path, for several temperatures. The points
are computed by taking the quotient between the spectral weight of the magnon
resonance, obtained from the $1/\mathcal{N}$-corrected
$S^{zz}(\bm{q},\omega)$, and the total (saddle-point plus
$1/\mathcal{N}$-corrections) spectral weight. We identify the crossover
temperature $T^{*}$ as the one above which, for all momenta except in the
vicinity of the Goldstone modes, the magnon spectral weight goes down to
10$\%$ of the total spectral weight. A similar behavior is observed in the
excitation spectrum of the superfluid 4He, where phonon excitations terminate
at some momentum by decaying into two roton excitations [38, 55] . Though the
decay mechanism differs, we term this antiferromagnetic region, between
$T^{*}$ and $T_{N}$, the terminated Goldstone regime (see Fig. 1). This
observation contrasts with the expectation that the magnon would persist
throughout the entire Brillouin zone for non-frustrated antiferromagnets.
Finally, at the Néel temperature, the whole spectra become incoherent, even at
the Goldstone modes.

Figure 7: Relative magnon spectral weight as a function of the temperature for
some momenta k on the Brillouin zone path
$K=(\frac{2\pi}{3},\frac{2\pi}{\sqrt{3}},\pi)\longrightarrow
M=(\pi,\pi/\sqrt{3},\pi)$. On the x-axis, the component common to all points,
$k_{z}=\pi$, is omitted.

###  IV.3 Lifetime of the magnon excitations

In order to further quantify the effect of the thermal activated spinons on
the two-spinon bound states, Fig. 8 shows the temperature evolution of its
width for a momentum $\bm{k}$ close to the Goldstone mode. While at $T=0$ the
finite width is given by twice the artificial broadening $\eta=0.02$,
necessary to perform the analytical continuation (see Eq. 24), as the
temperature rises the width has a $T^{2}$ behavior (see inset). At low
temperature ($T<T^{*}$) the lifetime hardly changes with respect to its ground
state value whereas for $T>T^{*}$ there is a pronounced increase. This is
consistent with the terminated Goldstone regime picture. It is worth to stress
that the thermal evolution of the magnon lifetime has been little explored in
the literature [56]. We recall that the present result is limited to two-
spinon scattering processes, and if we would like to include magnon-magnon
interactions, four spinon processes should be considered [30].

Figure 8: Temperature dependence of the quasiparticle peak width for momentum
$\textbf{k}=(2.225,3.4,\pi)$ which is close the Goldstone mode. Inset: Low
temperature fit of the form $f(T)=aT^{2}+b$, where $a=1.1555$ and $b=0.0408$.
$b$ is twice the artificial broadening $\eta$.

##  V Concluding remarks

Our results on the quasi-2D triangular antiferromagnet can be related with
recent numerical predictions for the specific heat and entropy of the 2D
triangular Heisenberg antiferromagnet [35, 34, 36]. In these studies two
characteristics temperatures are found: a high temperature scale $T_{h}\sim
0.5J$ signalling the development of nearest-neighbor spin correlation and a
low temperature scale $T_{l}\sim 0.2J$ related to the thermal activation of
the roton excitations. We conjecture that the crossover temperature
$T^{*}\simeq 0.75T_{N}$ can be associated with $T_{l}:$ around $T\sim T^{*}$
the antiferromagnet enters the terminated Goldstone regime, where the high
density of thermally activated spinons around $M$ gives rise to the decay of
the magnons away from the Goldstone points. On the other hand, the high
temperature scale $T_{h}$ can be associated with the ”pathological”
temperature $T\simeq 0.43J$ of the Schwinger boson theory [50, 49], above
which the mean-field solution corresponds to an infinite temperature perfect
paramagnet, in which the spin correlations vanish exactly.

In summary, in order to make contact with the experimental conditions we have
studied the thermal evolution of the magnetic excitations of a quasi-2D
triangular antiferromagnet, by means of a parton theory solved as a
$1/\mathcal{N}$ expansion of Schwinger boson theory. Recently, this approach
has been successfully applied –at zero temperature– to frustrated
antiferromagnets close to a quantum critical point [28, 33]. We have found
that the magnon excitations, described in our theory as two-spinon bound
states, decay due to the proliferation of thermally activated spinons. The
decay is very pronounced in the temperature range $T^{*}<T<T_{N}$, which we
have named terminated Goldstone regime, because –above $T^{*}$– the two-spinon
bound states survive only in the vicinity of the Goldstone modes. In the
quantum critical regime, above the Néel temperature, the excitations would
correspond to quasi-free spinon states.

The coexistence of quantum critical features and magnon excitations, below
Néel temperature, has been observed in weakly 3D interacting chains [37] and
also in the triangular antiferromagnet CuCrO2 [57]. These results gives
support to the idea that the fractionalization of magnons –proposed for 2D
frustrated AF– near a quantum critical point [31, 39, 4, 8] can be safely
extended to more realistic quasi-2D triangular antiferromagnets.

##  VI Acknowledgments

We thank Cristian Batista for insightful discussions. We acknowledge financial
support by CONICET (Argentina) under Grant PIP No. 3220.

##  Appendix A Numerical inverse of the dynamical matrix

The dynamical matrix (9) can be written as

$\mathcal{G}_{sp}^{-1}(\textbf{k},i\omega)=i\omega\mathcal{S}+\mathcal{M}_{\textbf{k}}$
(27)

where we have separated the k dependent terms from the $\omega$ dependent
ones, and $\mathcal{S}=diag(+1,-1,+1,-1)$ is the paraunitary matrix. Then we
perform a paraunitary diagonalization of the momentum dependent part,
$\mathcal{M}_{\textbf{k}}=\mathcal{U}^{\dagger}_{\textbf{k}}\mathcal{D}_{\textbf{k}}\mathcal{U}_{\textbf{k}}$,
where $\mathcal{D}_{\textbf{k}}$ is the diagonal matrix that contains the
eigenvalues and $\mathcal{U}_{\textbf{k}}$ is the paraunitary transformation
matrix [58]. Given that $\mathcal{S}$ is invariant under paraunitary
transformations
($S=\mathcal{U}^{\dagger}_{\textbf{k}}S\mathcal{U}_{\textbf{k}}$) we can write

$\mathcal{G}_{sp}^{-1}(\textbf{k},i\omega)=\mathcal{U}^{\dagger}_{\textbf{k}}\left(i\omega\mathcal{S}+\mathcal{D}_{\textbf{k}}\right)\mathcal{U}_{\textbf{k}}$
(28)

and the dynamical matrix inverse is obtained as

$\mathcal{G}_{sp}(\textbf{k},i\omega)=\mathcal{U}^{-1}_{\textbf{k}}\mathcal{S}\left(i\omega+\mathcal{S}\mathcal{D}_{\textbf{k}}\right)^{-1}\left(\mathcal{U}^{\dagger}_{\textbf{k}}\right)^{-1}$
(29)

where

$\left(i\omega+\mathcal{S}\mathcal{D}_{\textbf{k}}\right)^{-1}=\begin{pmatrix}\frac{1}{i\omega+\varepsilon^{+}_{\textbf{k}}}&0&0&0\\\
0&\frac{1}{i\omega-\varepsilon^{-}_{\textbf{k}}}&0&0\\\
0&0&\frac{1}{i\omega+\varepsilon^{-}_{\textbf{k}}}&0\\\
0&0&0&\frac{1}{i\omega-\varepsilon^{+}_{\textbf{k}}}\end{pmatrix}.$ (30)

The eigenvalues $\varepsilon^{\sigma}_{\textbf{k}}$ structure is consistent
with the $\textbf{k}\rightarrow-\textbf{k}$ invariance. Finally, the simple
fraction decomposition of the SP Green function matrix is

$\mathcal{G}_{sp}(\textbf{k},i\omega)=\sum_{\sigma\sigma^{\prime}}\frac{g^{\sigma\sigma^{\prime}}(\textbf{k})}{i\omega+\sigma\varepsilon^{\sigma^{\prime}}_{\textbf{k}}}$
(31)

where $\sigma,\sigma^{\prime}=\pm$ and

$\begin{split}g^{++}_{\alpha\beta}=\mathcal{U}^{-1}_{\alpha
1}\left(\mathcal{U}^{\dagger}\right)_{1\beta}^{-1},\hskip
14.22636ptg^{--}_{\alpha\beta}=-\mathcal{U}^{-1}_{\alpha
2}\left(\mathcal{U}^{\dagger}\right)_{2\beta}^{-1},\\\
g^{+-}_{\alpha\beta}=\mathcal{U}^{-1}_{\alpha
3}\left(\mathcal{U}^{\dagger}\right)_{3\beta}^{-1},\hskip
14.22636ptg^{-+}_{\alpha\beta}=-\mathcal{U}^{-1}_{\alpha
4}\left(\mathcal{U}^{\dagger}\right)_{4\beta}^{-1}.\end{split}$

In the last expressions we have omitted the k dependence to simplify the
notation.

## References

  * Laughlin and Pines [2000] R. B. Laughlin and D. Pines, The Theory of Everything, Proceedings of the National Academy of Sciences 97, 28 (2000).
  * Wen [2019] X.-G. Wen, Choreographed entanglement dances: Topological states of quantum matter, Science 363, 3099 (2019).
  * Anderson [1973] P. Anderson, Resonating valence bonds: A new kind of insulator?, Mater. Res. Bull. 8, 153 (1973).
  * Sachdev [2008] S. Sachdev, Quantum magnetism and criticality, Nature Physics 4, 173 (2008).
  * Normand [2009] B. Normand, Frontiers in frustrated magnetism, Contemporary Physics 50, 533 (2009).
  * Powell and McKenzie [2011] B. J. Powell and R. H. McKenzie, Quantum frustration in organic Mott insulators: from spin liquids to unconventional superconductors, Reports on Progress in Physics 74, 056501 (2011).
  * Savary and Balents [2017] L. Savary and L. Balents, Quantum spin liquids: a review, Rep. Prog. Phys. 80, 016502 (2017).
  * Zhou  _et al._ [2017] Y. Zhou, K. Kanoda, and T.-K. Ng, Quantum spin liquid states, Rev. Mod. Phys. 89, 025003 (2017).
  * Broholm  _et al._ [2020] C. Broholm, R. J. Cava, S. A. Kivelson, D. G. Nocera, M. R. Norman, and T. Senthil, Quantum spin liquids, Science 367, 0668 (2020).
  * Lee  _et al._ [2006] P. A. Lee, N. Nagaosa, and X.-G. Wen, Doping a Mott insulator: Physics of high-temperature superconductivity, Rev. Mod. Phys. 78, 17 (2006).
  * Kalmeyer and Laughlin [1987] V. Kalmeyer and R. B. Laughlin, Equivalence of the resonating-valence-bond and fractional quantum Hall states, Phys. Rev. Lett. 59, 2095 (1987).
  * Lacroix  _et al._ [2011] C. Lacroix, P. Mendels, and F. Mila, eds.,  _Introduction to Frustrated Magnetism: Materials, Experiments, Theory_ (Springer-Verlag, Berlin Heidelberg, 2011).
  * Kim  _et al._ [2019] T. Kim, K. Park, J. C. Leiner, and J.-G. Park, Hybridization and Decay of Magnetic Excitations in Two-Dimensional Triangular Lattice Antiferromagnets, Journal of the Physical Society of Japan 88, 081003 (2019).
  * Tennant [2019] A. Tennant, Studies of Spinons, Majoranas, and Monopoles in Spin Liquid and Quantum Critical Magnets with Neutrons, Journal of the Physical Society of Japan 88, 081009 (2019).
  * Tokura  _et al._ [2017] Y. Tokura, M. Kawasaki, and N. Nagaosa, Emergent functions of quantum materials, Nature Physics 13, 1056 (2017).
  * Keimer and Moore [2017] B. Keimer and J. E. Moore, The physics of quantum materials, Nature Physics 13, 1045 (2017).
  * Li [2019] Y. Li, YbMgGaO4: A Triangular-Lattice Quantum Spin Liquid Candidate, Advanced Quantum Technologies 2, 1900089 (2019).
  * Wen  _et al._ [2019] J. Wen, S.-L. Yu, S. Li, W. Yu, and J.-X. Li, Experimental identification of quantum spin liquids, npj Quantum Materials 4, 12 (2019).
  * Li  _et al._ [2020] Y. Li, P. Gegenwart, and A. A. Tsirlin, Spin liquids in geometrically perfect triangular antiferromagnets, Journal of Physics: Condensed Matter 32, 224004 (2020).
  * Doi  _et al._ [2004] Y. Doi, Y. Hinatsu, and K. Ohoyama, Structural and magnetic properties of pseudo-two-dimensional triangular antiferromagnets Ba3 MSb2O9 (M = Mn, Co, and Ni), J. Phys.: Condens. Matter 16, 8923 (2004).
  * Susuki  _et al._ [2013] T. Susuki, N. Kurita, T. Tanaka, H. Nojiri, A. Matsuo, K. Kindo, and H. Tanaka, Magnetization Process and Collective Excitations in the $S\mathbf{=}1/2$ Triangular-Lattice Heisenberg Antiferromagnet Ba3CoSb2O9, Phys. Rev. Lett. 110, 267201 (2013).
  * Koutroulakis  _et al._ [2015] G. Koutroulakis, T. Zhou, Y. Kamiya, J. D. Thompson, H. D. Zhou, C. D. Batista, and S. E. Brown, Quantum phase diagram of the ${S}=\frac{1}{2}$ triangular-lattice antiferromagnet Ba3CoSb2O9, Phys. Rev. B 91, 024410 (2015).
  * Ma  _et al._ [2016] J. Ma, Y. Kamiya, T. Hong, H. B. Cao, G. Ehlers, W. Tian, C. D. Batista, Z. L. Dun, H. D. Zhou, and M. Matsuda, Static and dynamical properties of the spin-$1/2$ equilateral triangular-lattice antiferromagnet Ba3CoSb2O9, Phys. Rev. Lett. 116, 087201 (2016).
  * Kamiya  _et al._ [2018] Y. Kamiya, L. Ge, T. Hong, Y. Qiu, D. Quintero-Castro, Z. Lu, H. Cao, M. Matsuda, E. S. Choi, C. Batista, M. Mourigal, H. D. Zhou, and J. Ma, The nature of spin excitations in the one-third magnetization plateau phase of Ba3CoSb2O9, Nature Communications 9 (2018).
  * Ito  _et al._ [2017] S. Ito, N. Kurita, H. Tanaka, S. Ohira-Kawamura, K. Nakajima, S. Itoh, K. Kuwahara, and K. Kakurai, Structure of the magnetic excitations in the spin-1/2 triangular-lattice Heisenberg antiferromagnet Ba3CoSb2O9, Nat. Comm. 8, 235 (2017).
  * Macdougal  _et al._ [2020] D. Macdougal, S. Williams, D. Prabhakaran, R. I. Bewley, D. J. Voneshen, and R. Coldea, Avoided quasiparticle decay and enhanced excitation continuum in the spin-$\frac{1}{2}$ near-Heisenberg triangular antiferromagnet ${\mathrm{ba}}_{3}{\mathrm{cosb}}_{2}{\mathrm{o}}_{9}$, Phys. Rev. B 102, 064421 (2020).
  * Chi  _et al._ [2022] R. Chi, Y. Liu, Y. Wan, H.-J. Liao, and T. Xiang, Spin excitation spectra of anisotropic spin-$1/2$ triangular lattice Heisenberg antiferromagnets, Phys. Rev. Lett. 129, 227201 (2022).
  * Ghioldi  _et al._ [2022] E. A. Ghioldi, S.-S. Zhang, Y. Kamiya, L. O. Manuel, A. E. Trumper, and C. D. Batista, Evidence of two-spinon bound states in the magnetic spectrum of Ba3CoSb2O9, Phys. Rev. B 106, 064418 (2022).
  * Ghioldi  _et al._ [2015] E. A. Ghioldi, A. Mezio, L. O. Manuel, R. R. P. Singh, J. Oitmaa, and A. E. Trumper, Magnons and excitation continuum in xxz triangular antiferromagnetic model: Application to Ba3CoSb2O9, Phys. Rev. B 91, 134423 (2015).
  * Ghioldi  _et al._ [2018] E. A. Ghioldi, M. G. Gonzalez, S.-S. Zhang, Y. Kamiya, L. O. Manuel, A. E. Trumper, and C. D. Batista, Dynamical structure factor of the triangular antiferromagnet: Schwinger boson theory beyond mean field, Phys. Rev. B 98, 184403 (2018).
  * Chubukov  _et al._ [1994] A. V. Chubukov, S. Sachdev, and T. Senthil, Quantum phase transitions in frustrated quantum antiferromagnets, Nuclear Physics B 426, 601 (1994).
  * Coldea  _et al._ [2003] R. Coldea, D. A. Tennant, and Z. Tylczynski, Extended scattering continua characteristic of spin fractionalization in the two-dimensional frustrated quantum magnet ${\mathrm{Cs}}_{2}{\mathrm{CuCl}}_{4}$ observed by neutron scattering, Phys. Rev. B 68, 134424 (2003).
  * Scheie  _et al._ [2024] A. O. Scheie, E. A. Ghioldi, J. Xing, J. A. M. Paddison, N. E. Sherman, M. Dupont, L. D. Sanjeewa, S. Lee, A. J. Woods, D. Abernathy, D. M. Pajerowski, T. J. Williams, S.-S. Zhang, L. O. Manuel, A. E. Trumper, C. D. Pemmaraju, A. S. Sefat, D. S. Parker, T. P. Devereaux, R. Movshovich, J. E. Moore, C. D. Batista, and D. A. Tennant, Proximate spin liquid and fractionalization in the triangular antiferromagnet KYbSe2, Nature Physics 20, 74 (2024).
  * Chen  _et al._ [2019] L. Chen, D.-W. Qu, H. Li, B.-B. Chen, S.-S. Gong, J. von Delft, A. Weichselbaum, and W. Li, Two-temperature scales in the triangular-lattice Heisenberg antiferromagnet, Phys. Rev. B 99, 140404 (2019).
  * Prelovsek and Kokalj [2018] P. Prelovsek and J. Kokalj, Finite-temperature properties of the extended Heisenberg model on a triangular lattice, Phys. Rev. B 98, 035107 (2018).
  * Gonzalez  _et al._ [2022] M. G. Gonzalez, B. Bernu, L. Pierre, and L. Messio, Ground-state and thermodynamic properties of the spin-$\frac{1}{2}$ Heisenberg model on the anisotropic triangular lattice, SciPost Phys. 12, 112 (2022).
  * Lake  _et al._ [2005] B. Lake, D. A. Tennant, C. D. Frost, and S. E. Nagler, Quantum criticality and universal scaling of a quantum antiferromagnet, Nature Materials 4, 329 (2005).
  * Stone  _et al._ [2006] M. B. Stone, I. A. Zaliznyak, T. Hong, C. L. Broholm, and D. H. Reich, Quasiparticle breakdown in a quantum spin liquid, Nature 440, 187 (2006).
  * Chubukov and Starykh [1995] A. V. Chubukov and O. A. Starykh, Confinement of spinons in the CPM-1 model, Phys. Rev. B 52, 440 (1995).
  * Arovas and Auerbach [1988] D. P. Arovas and A. Auerbach, Functional integral theories of low-dimensional quantum Heisenberg models, Phys. Rev. B 38, 316 (1988).
  * Auerbach [1994] A. Auerbach,  _Interacting electrons and quantum magnetism_ (Springer-Verlag, New York, 1994).
  * Trumper  _et al._ [1997] A. E. Trumper, L. O. Manuel, C. J. Gazza, and H. A. Ceccatto, Schwinger-Boson Approach to Quantum Spin Systems: Gaussian Fluctuations in the “Natural” Gauge, Phys. Rev. Lett. 78, 2216 (1997).
  * Zhang  _et al._ [2019] S.-S. Zhang, E. A. Ghioldi, Y. Kamiya, L. O. Manuel, A. E. Trumper, and C. D. Batista, Large-$s$ limit of the large-$n$ theory for the triangular antiferromagnet, Phys. Rev. B 100, 104431 (2019).
  * Ceccatto  _et al._ [1993] H. A. Ceccatto, C. J. Gazza, and A. E. Trumper, Nonclassical disordered phase in the strong quantum limit of frustrated antiferromagnets, Phys. Rev. B 47, 12329 (1993).
  * Flint and Coleman [2009] R. Flint and P. Coleman, Symplectic $N$ and time reversal in frustrated magnetism, Phys. Rev. B 79, 014424 (2009).
  * Mezio  _et al._ [2011] A. Mezio, C. N. Sposetti, L. O. Manuel, and A. E. Trumper, A test of the bosonic spinon theory for the triangular antiferromagnet spectrum, EPL (Europhysics Letters) 94, 47001 (2011).
  * Azaria  _et al._ [1993] P. Azaria, B. Delamotte, and D. Mouhanna, Spontaneous symmetry breaking in quantum frustrated antiferromagnets, Phys. Rev. Lett. 70, 2483 (1993).
  * Yoshioka and Miyazaki [1991] D. Yoshioka and J. Miyazaki, Boson Mean Field Theory of the Triangular Lattice Heisenberg Model, Journal of the Physical Society of Japan 60, 614 (1991).
  * De Silva  _et al._ [2002] T. N. De Silva, M. Ma, and F. C. Zhang, Pathology of Schwinger boson mean-field theory for Heisenberg spin models, Phys. Rev. B 66, 104417 (2002).
  * Tchernyshyov and Sondhi [2002] O. Tchernyshyov and S. L. Sondhi, Liquid gas and other unusual thermal phase transitions in some large-$N$ magnets, Nucl. Phys. B 639, 429 (2002), arXiv:cond-mat/0202128 .
  * A. Du [2002] J. L. A. Du, G.Z. Wei, Low-temperature properties of the Heisenberg antiferromagnets on the layered triangular lattice, phys. stat. sol. (b) 234, 636 (2002).
  * Yasuda  _et al._ [2005] C. Yasuda, S. Todo, K. Hukushima, F. Alet, M. Keller, M. Troyer, and H. Takayama, Néel Temperature of Quasi-Low-Dimensional Heisenberg Antiferromagnets, Phys. Rev. Lett. 94, 217201 (2005).
  * Manuel  _et al._ [1998] L. O. Manuel, A. E. Trumper, and H. A. Ceccatto, Rotational invariance and order-parameter stiffness in frustrated quantum spin systems, Phys. Rev. B 57, 8348 (1998).
  * Mezio  _et al._ [2012] A. Mezio, L. O. Manuel, R. R. P. Singh, and A. E. Trumper, Low temperature properties of the triangular-lattice antiferromagnet: a bosonic spinon theory, New Journal of Physics 14, 123033 (2012).
  * Zhitomirsky and Chernyshev [2013] M. E. Zhitomirsky and A. L. Chernyshev, Colloquium : Spontaneous magnon decays, Rev. Mod. Phys. 85, 219 (2013).
  * Bayrakci  _et al._ [2013] S. P. Bayrakci, D. A. Tennant, P. Leininger, T. Keller, M. C. R. Gibson, S. D. Wilson, R. J. Birgeneau, and B. Keimer, Lifetimes of antiferromagnetic magnons in two and three dimensions: Experiment, theory, and numerics, Phys. Rev. Lett. 111, 017204 (2013).
  * Kajimoto  _et al._ [2015] R. Kajimoto, K. Tomiyasu, K. Nakajima, S. Ohira-Kawamura, Y. Inamura, and T. Okuda, Development of Spin Correlations in the Geometrically Frustrated Triangular-Lattice Heisenberg Antiferromagnet CuCrO2, Journal of the Physical Society of Japan 84, 074708 (2015).
  * Colpa [1978] J. H. P. Colpa, Diagonalization of the quadratic boson Hamiltonian, Physica A 93, 327 (1978).