# Thermal squeezing and nonlinear spectral shift of magnons in
antiferromagnetic insulators

Mahroo Shiranzaei  Division of Materials Theory, Department of Physics and
Astronomy, Uppsala University, Box 516, SE-75120 Uppsala, Sweden  Roberto E.
Troncoso  Center for Quantum Spintronics, Department of Physics, Norwegian
University of Science and Technology, NO-7491 Trondheim, Norway  Department of
Mechanical and Industrial Engineering, Norwegian University of Science and
Technology, Høgskoleringen 1, 7491 Trondheim, Norway  Jonas Fransson  Division
of Materials Theory, Department of Physics and Astronomy, Uppsala University,
Box 516, SE-75120 Uppsala, Sweden  Arne Brataas  Center for Quantum
Spintronics, Department of Physics, Norwegian University of Science and
Technology, NO-7491 Trondheim, Norway  Alireza Qaiumzadeh  Center for Quantum
Spintronics, Department of Physics, Norwegian University of Science and
Technology, NO-7491 Trondheim, Norway

###### Abstract

We investigate the effect of magnon-magnon interactions on the dispersion and
polarization of magnons in collinear antiferromagnetic (AF) insulators at
finite temperatures. In two-sublattice AF systems with either uniaxial or
biaxial magnetocrystalline anisotropies, we implement a self-consistent
Hartree-Fock mean-field approximation to explore the nonlinear thermal
interactions. The resulting nonlinear magnon interactions separate into two-
magnon intra- and interband scattering processes. Furthermore, we compute the
temperature dependence of the magnon bandgap and AF resonance modes due to
nonlinear magnon interactions for square and hexagonal lattices. In addition,
we study the effect of magnon interactions on the polarization of magnon
modes. We find that although the noninteracting eigenmodes in the uniaxial
case are circularly polarized, but in the presence of nonlinear thermal
interactions the local U(1) symmetry of the Hamiltonian is broken. The
attractive nonlinear interactions squeeze the low energy magnon modes and make
them elliptical. In the biaxial case, on the other hand, the bare eigenmodes
of low energy magnons are elliptically polarized but thermal nonlinear
interactions squeeze them further. Direct measurements of the predicted
temperature-dependent AF resonance modes and their polarization can be used as
a tool to probe the nonlinear interactions. Our findings establish a framework
for exploring the effect of thermal magnon interactions in technologically
important magnetic systems, such as magnetic stability of recently discovered
two-dimensional magnetic materials, coherent transport of magnons, Bose-
Einstein condensation of magnons, and magnonic topological insulators.

##  I Introduction

Antiferromagnetic (AF) systems occupy a central position at the frontiers of
modern spintronics due to their functional potential in future energy-
efficient and ultrafast information and communication technologies Marti  _et
al._ (2014); Kriegner  _et al._ (2016); Lebrun  _et al._ (2018); Němec  _et
al._ (2018); Geprags  _et al._ (2020). AF ordered systems are abundant
materials compared to ferromagnets and range from collinear to noncollinear
and uniaxial to multi-axis materials. Current research centers on exploiting
their most remarkable properties as the lack of parasitic stray fields, THz
spin dynamics, and two-magnon-degrees of freedom Baltz  _et al._ (2018);
Qaiumzadeh  _et al._ (2018); Chumak  _et al._ (2014); Bowlan  _et al._ (2018);
Jungwirth  _et al._ (2016); Proskurin  _et al._ (2017). Of particular interest
are the electrical insulator class of AF materials Chumak  _et al._ (2015);
Gomonay and Loktev (2014); Daniels  _et al._ (2015); Khymyn  _et al._ (2016);
Jungwirth  _et al._ (2016); Gomonay  _et al._ (2017); Nishitani  _et al._
(2010, 2012); Lebrun  _et al._ (2020), where Joule-heating losses caused by
itinerant electrons are absent. Magnons, the low-energy bosonic excitations of
magnetically ordered materials, can carry spin angular momentum in AF
insulators across large distances, as recently shownLebrun  _et al._ (2020,
2018); Rezende  _et al._ (2019).

Magnonics is an emerging field in spintronics and magnetic-based
nanotechnology Chumak  _et al._ (2015) that exhibits various features from
quantum coherent phenomena Bunkov and Volovik (2010); Dzyapko  _et al._
(2017); Tupitsyn  _et al._ (2008); Rezende (2009), chiral magnonic spin
transport, magnonic crystals Troncoso  _et al._ (2015); Zakeri (2020); Chumak
_et al._ (2017); Singh  _et al._ (2021), and ultrafast magnetization dynamics
Kirilyuk  _et al._ (2010). Nonlinear magnon interactions are ubiquitous and
inevitable in magnetically ordered materials. Although in classical magnets
nonlinear magnon-magnon interactions are weak at low temperatures, their
effects significantly influence various phenomena. These phenomena are
relevant in spin transport Wang  _et al._ (2020); Bayrakci  _et al._ (2013)
through magnon relaxation and spin conductivity, as well as close to magnetic
phase transitions. Moreover, nonlinear magnon interactions are essential to
stabilize magnetic droplets Mohseni  _et al._ (2013) and Bose-Einstein
condensation of magnons Mohseni  _et al._ (2020); Tupitsyn  _et al._ (2008);
Arakawa (2019); Giamarchi  _et al._ (2008). It is also an important ingredient
in emergent magnonic topological phases Mook  _et al._ (2021). A recent study
demonstrated that nonlinear magnon interactions of Dirac magnons in honeycomb
ferromagnets strongly renormalize the magnon band structure and magnon
lifetimes Pershoguba  _et al._ (2018). In AF systems, nonlinear magnon
interactions also have significant effects even at low temperatures due to
zero-point fluctuations Zhitomirsky and Chernyshev (2013). Additionally, two-
magnon spectra of AF systems are strongly influenced by magnon-magnon
interactions Elliott and Thorpe (1969).

Nonlinear magnon interactions, even in simple AF systems, are pretty
complicated. Although there are many theoretical studies on this topic, e.g.,
see a recent review in Ref. Rezende  _et al._ , 2019, to the best of our
knowledge, a systematic investigation of nonlinear magnon interactions at
finite temperature in biaxial AF materials remains absent Liu (1992). On the
other hand, in AF systems, the underlying crystal lattice structure affects
the magnetic ground states and interactions. Therefore, one needs to study
each case separately Chernyshev and Zhitomirsky (2009).

The recent discovery of two-dimensional (2D) materials establishes a new
platform for the study of nonlinear magnon interactions in different lattice
structures Bayrakci  _et al._ (2013); McCreary  _et al._ (2020); Niyazi  _et
al._ (2021); Mkhitaryan and Ke (2021); Li  _et al._ (2018). In low-dimensional
magnetic systems, the presence of a finite bandgap in the low-energy magnon
dispersion is an essential ingredient for the stability of long-range magnetic
order at finite temperature and overcoming the bottleneck imposed by the
Mermin-Wagner-Hohenberg theorem Mermin and Wagner (1966); Hohenberg (1967).
Therefore, the effect of nonlinear thermal magnon interactions on the magnetic
resonance frequency and magnon dispersion is a key issue in the study of long-
range order stability of 2D magnetic systems.

There are already several primary works which investigate higher order magnon
interactions in AF systems within a mean-field approach. In Refs. Castilla and
Chakravarty, 1991; Igarashi, 1992, the effect of nonlinear interactions on
magnon dispersion of a square lattice Heisenberg antiferromagnet were studied
via a Green’s function method at zero temperature. In Ref. An  _et al._ , 2001
effect of nonlinear thermal magnon interactions on spin-wave excitation and
sublattice magnetization in a 3D anisotropic quantum Heisenberg
antiferromagnet with a tuneable interplanar coupling has been studied via an
equation of motion formalism within a mean-field approach.

In this paper, we first develop a formalism to compute the nonlinear
interactions of thermal magnons in both uniaxial and biaxial AF insulators in
arbitrary dimensions using a self-consistent Hartree-Fock mean-field theory.
Next, we investigate the effect of thermal four-magnon interactions on bare
magnon dispersion and AF resonance frequencies, which are directly accessible
experimentally Gitgeatpong  _et al._ (2017). Later, we investigate the effect
of nonlinear interactions on the polarization of magnon modes. For numerical
calculations, We consider collinear AF materials with uniaxial and biaxial
magnetic anisotropies in square and hexagonal lattices. We study these systems
as prototypes of recently discovered 2D AF materials to explore the stability
of these systems against thermal fluctuations.

We have organized this paper as follows. In Sec. II, we introduce the spin
Hamiltonian and review the linear spin-wave theory. In Sec. III, we present a
Hartree-Fock mean-field theory for thermal magnon interactions and compute
mean-field coefficients. In Sec. IV, we compute the nonlinear spectral shift
of magnons in uniaxial and biaxial AF systems and present numerical
calculations for two specific lattice models. In Sec. V, we investigate the
effect of nonlinear thermal interactions on the polarization of magnons. We
conclude in Sec. VI with a discussion of our results. In the appendices, we
present the technical details of the calculations.

##  II Spin Hamiltonian and its bosonic representation

We describe a generic two-sublattice AF insulator with the local spin field
$\mathbf{S}_{\mathcal{A}(\mathcal{B})}$ in sublattice $\mathcal{A}$
($\mathcal{B}$). The effective AF spin Hamiltonian reads, Baltz  _et al._
(2018); Rezende  _et al._ (2019),

$\displaystyle H_{s}=J\sum_{\langle
i,j\rangle}\mathbf{S}_{i}\cdot\mathbf{S}_{j}+\sum_{i}\mathcal{K}_{x}(S^{x}_{i})^{2}-\sum_{i}\mathcal{K}_{z}(S^{z}_{i})^{2},$
(1)

where $\langle i,j\rangle$ represents summation over all nearest neighbor
sites and the local spins $\mathbf{S}_{i\in\\{\mathcal{A},\mathcal{B}\\}}$
have the same amplitude $S$ at each lattice site. The Hamiltonian includes the
isotropic quantum AF Heisenberg exchange interaction, $J>0$, and single-ion-
type magnetic anisotropies, characterized by $\mathcal{K}_{x}>0$ and
$\mathcal{K}_{z}>0$, which represent the hard and easy magnetic axes,
respectively. In the uniaxial case with an easy axis, i.e.,
$\mathcal{K}_{x}=0$, both the spin Hamiltonian and the ground state have
$U(1)$ symmetry while in the biaxial case this symmetry is broken. The absence
of $U(1)$ symmetry in the ground state breaks the degeneracy of the two AF
magnon eigenmodes. We emphasize that although the spin Hamiltonian, Eq. (1),
has $U(1)$ symmetry even in the easy-plane case, i.e., $\mathcal{K}_{z}=0$,
this symmetry is spontaneously broken in the ground state. Thus, the
eigenmodes of the easy-axis AF systems are not degenerate in the later case
Qaiumzadeh  _et al._ (2017).

To study spin excitations, it is more convenient to rewrite the spin
Hamiltonian, Eq. (1), in terms of so-called deviation bosonic operators using
the Holstein-Primakoff (HP) transformation Holstein and Primakoff (1940). We
assume the ground state of the AF spin Hamiltonian is a collinear Néel state
with ${\bf S}_{i\in\\{{\cal A}\\}}=(0,0,S)$ and ${\bf S}_{j\in\\{{\cal
B}\\}}=(0,0,-S)$; and thus the HP transformation in the two-sublattice AF
system reads,

$\displaystyle S^{z}_{i}=S-a^{\dagger}_{i}a_{i},\hskip
7.11317ptS^{-}_{i}=a^{\dagger}_{i}\sqrt{2S-a^{\dagger}_{i}a_{i}},$ (2a)
$\displaystyle S^{z}_{j}=-S+b^{\dagger}_{j}b_{j},\hskip
2.84544ptS^{-}_{j}=\sqrt{2S-b^{\dagger}_{j}b_{j}}\;b_{j},$ (2b)

and $S^{+}_{i(j)}=\big{(}S^{-}_{i(j)}\big{)}^{\dagger}$, where $i$ and $j$
denote sites in the sublattices $\mathcal{A}$ and $\mathcal{B}$, respectively.
The deviation operators $a$ ($b$) and $a^{\dagger}$ ($b^{\dagger}$) are,
respectively, annihilation and creation operators at the sublattice
$\mathcal{A}$ ($\mathcal{B}$) that follow the bosonic communication relations.
The number of HP bosons at each site on sublattices $\mathcal{A}$ and
$\mathcal{B}$ is upper-limited by $2S$ König and Hucht (2021). Inserting Eqs.
(2a) and (2b) into the AF spin Hamiltonian, Eq. (1), and under the assumption
of small deviation number $\langle a^{\dagger}_{i}a_{i}\rangle,\langle
b^{\dagger}_{i}b_{i}\rangle\ll 2S$, with $\langle\cdot\rangle$ the expectation
value, we do a Taylor series expansion in the powers of $1/S$. This procedure
leads to an effective bosonic Hamiltonian with an arbitrary interaction order,
$H_{b}=E^{c}_{0}+H_{b}^{(2)}+H_{b}^{(4)}+\cdots$, where ellipsis stands for
higher order interactions. The classical ground-state energy of Eq. (1) is
$E^{c}_{0}=NS\left[-(zJ/2+\mathcal{K}_{z})S+\mathcal{K}_{x}/2)\right]$, where
$z$ denotes the coordination number and $N$ is the total number of sites, and
$H_{b}^{(2)}$ and $H_{b}^{(4)}$ represent the bosonic excitations.
$H_{b}^{(2)}$ consists of quadratic terms and describe the dynamics of
noninteracting bosonic quasiparticles while $H_{b}^{(4)}$ consists of quartic
terms and describes the two quasiparticle scatterings via lowest-order
nonlinear interactions, i.e., four-boson interactions. In this model, Eq. (1),
three-magnon interactions are forbidden and we do not have cubic terms
Zhitomirsky and Chernyshev (2013).

###  II.1 Linear spin-wave theory

At a low density of magnon excitations, we can disregard higher order
nonlinear magnon interactions. The noninteracting bosonic Hamiltonian
$H_{b}^{(2)}$ becomes diagonal in the momentum space after a canonical
Bogoliubov transformation by introducing new bosonic operators
$\alpha_{\boldsymbol{q}}$ and $\beta_{\boldsymbol{q}}$ at wavevector
${\boldsymbol{q}}$, see Appendix A for details. The resulting magnon
Hamiltonian in the eigenbasis is $H_{m}=E_{0}+{H}_{m}^{(2)}$. The ground-state
energy is $E_{0}=E^{c}_{0}+E^{q}_{0}$, where $E^{q}_{0}$ denotes quantum
correction to the classical ground-state energy, see Appendix A.2. The quantum
zero-point energy $E^{q}_{0}$ lowers the classical ground-state energy and
expresses the well-known fact that the true ground state of an AF system is
not a Néel state Auerbach (1998). The Hamiltonian of noninteracting magnons is
given by

$\displaystyle{H}_{m}^{(2)}=\sum_{\boldsymbol{q}}\left[\epsilon^{\alpha}_{\boldsymbol{q}}\alpha^{\dagger}_{\boldsymbol{q}}\alpha_{\boldsymbol{q}}+\epsilon^{\beta}_{\boldsymbol{q}}\beta^{\dagger}_{\boldsymbol{q}}\beta_{\boldsymbol{q}}\right].$
(3)

The bare eigenenergy for the magnon mode $\chi\in\\{\alpha,\beta\\}$ is given
by

$\displaystyle\epsilon^{\chi}_{\boldsymbol{q}}=S\sqrt{\big{(}zJ+\mathcal{K}_{x}+2\mathcal{K}_{z}\big{)}^{2}-\big{(}zJ|\gamma_{\boldsymbol{q}}|\mp\mathcal{K}_{x}\big{)}^{2}},$
(4)

where the sign $-$($+$) refers to $\alpha$($\beta$) magnon mode.
$\gamma_{\boldsymbol{q}}=z^{-1}\sum_{i=1}^{z}e^{i{\boldsymbol{q}}\cdot\boldsymbol{\delta}_{i}}$
is the structure factor, $\boldsymbol{\delta}_{i}$ denotes the nearest-
neighbor vectors. In the uniaxial easy-axis limit $\mathcal{K}_{x}=0$, and in
the absence of external magnetic fields and Dzyaloshinskii-Moriya interactions
(DMIs), the two magnon eigenenergies are degenerate and the eigenmodes are
circularly polarized with opposite chiralities Rezende  _et al._ (2019);
Qaiumzadeh  _et al._ (2018). The presence of a hard-axis anisotropy
($\mathcal{K}_{x}\neq 0$) removes the degeneracy of magnon eigenstates and
causes the two chiral eigenmodes become elliptically polarized Rezende  _et
al._ (2019); Qaiumzadeh  _et al._ (2018). The latter is a consequence of the
breaking of $U(1)$ symmetry in the ground state of the Hamiltonian 1. The
structure factors for 2D square ($z=4$) and hexagonal ($z=3$) lattices are,
$\gamma_{\boldsymbol{q}}=2\left(\cos(q_{x}/q_{m})+\cos(q_{y}/q_{m})\right)/z$
and
$\gamma_{\boldsymbol{q}}=e^{iq_{x}/q_{m}}\left(1+2e^{-i3q_{x}/2q_{m}}\cos(\sqrt{3}q_{y}/2q_{m})\right)/z$,
respectively, where $q_{m}=1/a$ and $a$ is the lattice constant. In the long-
wavelength limit ${\bm{q}}\rightarrow 0$, the structure factors of both
hexagonal and square lattices reduce to $\gamma_{\boldsymbol{q}}\simeq
1-q^{2}/4q_{m}^{2}$; and thus the mode-dependent eigenenergies becomes to
$\epsilon^{\alpha}_{\mathbf{q}}=S\sqrt{4(\mathcal{K}_{z}+\mathcal{K}_{x})(zJ+\mathcal{K}_{z})+zJq^{2}(zJ-\mathcal{K}_{x})/2q_{m}^{2}}$
and
$\epsilon^{\beta}_{\mathbf{q}}=S\sqrt{4\mathcal{K}_{z}(zJ+\mathcal{K}_{z}+\mathcal{K}_{x})+zJq^{2}(zJ+\mathcal{K}_{x})/2q_{m}^{2}}$.
In other words, around the center of the Brillouin zone, called
$\Gamma$-point, the two low-energy AF magnon modes propagate as relativistic-
like quasiparticles with different effective _rest masses_ and _speed of
light_ , determined by the magnetic exchange stiffness and magnetic
anisotropies. The bare magnon bandgaps of each mode are
$\epsilon^{\alpha}_{\mathbf{q}=0}$ and $\epsilon^{\beta}_{\mathbf{q}=0}$ that
can be measured by in AF resonance experiments. The presence of this gap
stabilizes long-range magnetic order at finite temperature in 2D magnetic
systems and the magnetic transition temperatures in these systems are
proportional to this energy gap. We will later show how the bare magnon
dispersions become renormalized in the presence of nonlinear thermal magnon
interactions.

##  III Mean-Field theory of Magnon-Magnon Interactions

In this section, we use a Hartree-Fock mean-field theory to treat the
nonlinear thermal magnon interactions in a two-sublattice AF system at finite
temperature. The lowest-order nonlinear terms related to two-boson scattering
are represented by four-boson interactions,

$\displaystyle H^{(4)}_{b}=\frac{-J}{4}\sum_{\langle
i,j\rangle}\Big{[}a^{\dagger}_{i}a_{i}b^{\dagger}_{j}b_{j}+a_{i}b^{\dagger}_{j}b_{j}b_{j}+(a\leftrightarrow
b)\Big{]}-\frac{\mathcal{K}_{x}}{8}\sum_{i}\Big{[}\Delta\mathcal{K}a^{\dagger}_{i}a^{\dagger}_{i}a_{i}a_{i}+a^{\dagger}_{i}a_{i}a_{i}a_{i}+(a\leftrightarrow
b)\Big{]}+\text{H.c.},$ (5)

with $\Delta\mathcal{K}=2(1+\mathcal{K}_{z}/\mathcal{K}_{x})$, and H.c. is the
Hermitian conjugate. Four-boson interactions in AF systems are more
complicated than the ferromagnetic one. In ferromagnets, in the absence of
long-range dipolar fields, the four-boson interactions conserve the number of
quasiparticles because there are only quartic terms with two creation and two
annihilation operators. On the contrary, as it can be seen in Eq. (5), in AF
systems the magnon number could in general be a nonconserved quantity. It was
also shown that spontaneous three-magnon and in general all other $n$-magnon
decays are forbidden in collinear AF systems by energy conservation at zero
temperature in the absence of an external magnetic field Harris  _et al._
(1971); Zhitomirsky and Chernyshev (2013).

To explore the effects of nonlinear magnon interactions on the bare spectra of
magnons at finite temperature, we implement a mean-field theory based on the
self-consistent Hartree-Fock approximation Balucani  _et al._ (1980); Oguchi
(1960) to reduce Eq. (5) into an effective noninteracting Hamiltonian in terms
of quadratic magnon operators. In the mean-field treatment, we start out by
substituting
$\xi^{\dagger}_{i}{\xi}^{\prime}_{j}=\langle\xi^{\dagger}_{i}{\xi}^{\prime}_{j}\rangle_{\text{th}}+\Xi_{\xi{\xi}^{\prime}}$
into the Hamiltonian 5, where
$\xi_{i},{\xi}^{\prime}_{i}\in\\{a_{i},b_{i},a_{i}^{\dagger},b_{i}^{\dagger}\\}$
and $\langle\cdot\rangle_{\text{th}}$ is a thermal average. We expand up to
the first order in the field
$\Xi_{\xi{\xi}^{\prime}}\equiv\xi^{\dagger}_{i}{\xi}^{\prime}_{j}-\langle\xi^{\dagger}_{i}{\xi}^{\prime}_{j}\rangle_{\text{th}}$,
and express the operators in the basis of $\alpha$\- and $\beta$-magnons
through the Bogoliubov transformation, see the Appendix A for details. In the
leading-order correction, it is assumed that mean-field parameter
$\Xi_{\xi{\xi}^{\prime}}$ characterizes slight nonequilibrium deviations from
the thermal equilibrium, which is consistent with our earlier assumption of
small deviation in the boson number in the Taylor series expansion of HP
transformations. To guarantee the thermalization of each magnon mode, we
assume that nonlinear magnon scattering-times are much faster than the
timescales associated with spin non-conserving interactions such as magnon-
phonon interactions Troncoso  _et al._ (2020). In the Hartree-Fock
approximation, we only keep the diagonal terms and assume for each bosonic
operator
$\chi_{\boldsymbol{q}}\in\\{\alpha_{\boldsymbol{q}},\beta_{\boldsymbol{q}}\\}$,
so that
$\langle\chi^{\dagger}_{\boldsymbol{q}}\chi_{\boldsymbol{q}^{\prime}}\rangle_{\text{th}}=\delta_{{\boldsymbol{q}},{\boldsymbol{q}}^{\prime}}n_{\chi,{\boldsymbol{q}}}$
Li  _et al._ (2018); Mkhitaryan and Ke (2021). The thermodynamic equilibrium
of the AF magnon number at each magnon mode $\chi\in\\{\alpha,\beta\\}$ is
defined through the Bose-Einstein distribution function
$n_{\chi,{\boldsymbol{q}}}=(e^{\epsilon_{\chi,{\boldsymbol{q}}}/k_{B}T}-1)^{-1}$,
with the Boltzmann constant $k_{B}$, and the temperature $T$. The chemical
potential of magnons is zero at thermal equilibrium Flebus (2019).

The effective mean-field Hamiltonian in momentum space and first order in the
fluctuations $\Xi_{\xi{\xi}^{\prime}}$, becomes
${H}_{\text{MF}}^{(4)}=\frac{1}{2}\sum_{\boldsymbol{q}}\Psi^{\dagger}_{\boldsymbol{q}}\mathbb{H}^{(4)}\Psi_{\boldsymbol{q}}$,
with the vector operator
$\Psi_{\boldsymbol{q}}=(\alpha_{\boldsymbol{q}},\beta_{\boldsymbol{q}},\alpha^{\dagger}_{-\boldsymbol{q}},\beta^{\dagger}_{-\boldsymbol{q}})^{T}$.
The matrix Hamiltonian $\mathbb{H}^{(4)}$ is defined by,

$\displaystyle\mathbb{H}^{(4)}=\left(\begin{array}[]{cccc}{\cal
C}_{\alpha\alpha}&{\cal D}_{\alpha\beta}&2\bar{\cal
C}^{\prime}_{\alpha\alpha}&\bar{\cal D}^{\prime}_{\alpha\beta}\\\ \bar{\cal
D}_{\alpha\beta}&{\cal C}_{\beta\beta}&\bar{\cal
D}^{\prime}_{\alpha\beta}&2\bar{\cal C}^{\prime}_{\beta\beta}\\\ 2{\cal
C}^{\prime}_{\alpha\alpha}&{\cal D}^{\prime}_{\alpha\beta}&{\cal
C}_{\alpha\alpha}&\bar{\cal D}_{\alpha\beta}\\\ {\cal
D}^{\prime}_{\alpha\beta}&2{\cal C}^{\prime}_{\beta\beta}&{\cal
D}_{\alpha\beta}&{\cal C}_{\beta\beta}\end{array}\right),$ (10)

where the overbar represents the complex conjugate and the explicit dependence
of the temperature-dependent mean-field coefficients on the wavevector and
temperature is omitted for brevity. The Hamiltonian ${H}_{\text{MF}}^{(4)}$
represents the quantum and temperature-dependent corrections to the
noninteracting magnon Hamiltonian, when the lowest-order nonlinear AF dynamics
is taken into account. The coefficients of the matrix $\mathbb{H}^{(4)}$
depend on temperature, for more details see Appendix B.2. This effective
Hamiltonian ${H}_{\text{MF}}^{(4)}$ describes scattering processes between
magnons of the same band (intraband contribution) and magnons of different
bands (interband contribution). The components ${\cal C}_{\chi\chi}$ and
${\cal C}^{\prime}_{\chi\chi}$, describing intraband processes of magnons at
the $\chi$-band, are given by Eqs. (30) in App. (B.2).

$\displaystyle{\cal C}_{\alpha\alpha,{\boldsymbol{q}}}$
$\displaystyle=\Delta^{-}_{\boldsymbol{q}}|u_{\boldsymbol{q}}|^{2}+\Delta^{-}_{-\boldsymbol{q}}|x_{\boldsymbol{q}}|^{2}+2\,\text{Re}\left[u_{\boldsymbol{q}}x_{\boldsymbol{q}}\,\big{(}\xi-\text{Re}\,\Gamma_{q}\big{)}\right],$
(11a) $\displaystyle{\cal C}_{\beta\beta,{\boldsymbol{q}}}$
$\displaystyle=\Delta^{+}_{\boldsymbol{q}}|v_{\boldsymbol{q}}|^{2}+\Delta^{+}_{\boldsymbol{q}}|w_{\boldsymbol{q}}|^{2}+2\,\text{Re}\left[v_{\boldsymbol{q}}w_{\boldsymbol{q}}\,\big{(}\xi+\text{Re}\,\Gamma_{q}\big{)}\right],$
(11b) $\displaystyle{\cal C}^{\prime}_{\alpha\alpha,{\boldsymbol{q}}}$
$\displaystyle=2\Delta^{-}_{\boldsymbol{q}}u_{\boldsymbol{q}}{x}_{\boldsymbol{q}}+\Lambda^{-}_{\boldsymbol{q}}u^{2}_{\boldsymbol{q}}+\bar{\Lambda}^{-}_{\boldsymbol{q}}{x}^{2}_{\boldsymbol{q}},$
(11c) $\displaystyle{\cal C}^{\prime}_{\beta\beta,{\boldsymbol{q}}}$
$\displaystyle=2\Delta^{+}_{\boldsymbol{q}}v_{\boldsymbol{q}}{w}_{\boldsymbol{q}}+\Lambda^{+}_{\boldsymbol{q}}v^{2}_{\boldsymbol{q}}+\bar{\Lambda}^{+}_{\boldsymbol{q}}{w}^{2}_{\boldsymbol{q}},$
(11d)

where $\Delta^{\pm}_{\boldsymbol{q}}$, $\xi$, $\Gamma_{\boldsymbol{q}}$ and
$\Lambda^{\pm}_{\boldsymbol{q}}$ are defined in Eqs. (31). The Bogoliubov
coefficients $u_{\boldsymbol{q}},v_{\boldsymbol{q}},x_{\boldsymbol{q}}$ and
$w_{\boldsymbol{q}}$ are defined by Eqs. (22a)-(22d) in Appendix A. Interband
scattering of magnons are characterized by the components ${\cal
D}_{\alpha\beta}$ and ${\cal D}^{\prime}_{\alpha\beta}$, and read,

$\displaystyle{\cal D}_{\alpha\beta,{\boldsymbol{q}}}=$
$\displaystyle\Gamma_{\boldsymbol{q}}x_{\boldsymbol{q}}v_{\boldsymbol{q}}-\bar{\Gamma}_{\boldsymbol{q}}{u}_{\boldsymbol{q}}{w}_{\boldsymbol{q}}+\Theta_{-\boldsymbol{q}}x_{\boldsymbol{q}}{w}_{\boldsymbol{q}}-\Theta_{\boldsymbol{q}}v_{\boldsymbol{q}}{u}_{\boldsymbol{q}}$
$\displaystyle+2\,\xi
v_{\boldsymbol{q}}\,x_{\boldsymbol{q}}+2\,\bar{\xi}\,u_{\boldsymbol{q}}\,w_{\boldsymbol{q}},$
(12a) $\displaystyle{\cal D}^{\prime}_{\alpha\beta,{\boldsymbol{q}}}=$
$\displaystyle\bar{\Gamma}_{\boldsymbol{q}}{x}_{\boldsymbol{q}}{w}_{\boldsymbol{q}}-\Gamma_{\boldsymbol{q}}u_{\boldsymbol{q}}v_{\boldsymbol{q}}-{\Theta}_{-\boldsymbol{q}}{x}_{\boldsymbol{q}}v_{\boldsymbol{q}}+\Theta_{\boldsymbol{q}}{w}_{\boldsymbol{q}}u_{\boldsymbol{q}}$
$\displaystyle+2\xi^{\prime}u_{\boldsymbol{q}}v_{\boldsymbol{q}}+2\,\bar{\xi}^{\prime}\,x_{\boldsymbol{q}}w_{\boldsymbol{q}},$
(12b)

where $\Gamma^{\prime}_{\boldsymbol{q}}$, $\Theta_{\boldsymbol{q}}$ and
$\xi^{\prime}$ can be found in Eqs. (31). Through the interband process,
$\alpha$\- and $\beta$-magnon modes become coupled and thus, hybridized. The
definition of $A_{a},\,A^{\prime}_{a},\,A_{ab,{\boldsymbol{q}}}$ and
$A^{\prime}_{ab,{\boldsymbol{q}}}$ are respectively detailed in Eqs.
(29b)-(29e) in the Appendix B. The mean-field coefficients ${\cal
C}_{\chi\chi},\,{\cal C}^{\prime}_{\chi\chi},\,{\cal D}_{\alpha\beta}$ and
${\cal D}^{\prime}_{\alpha\beta}$ are functions of temperature and magnon
energy through the Bose-Einstein distribution $n_{\chi,\boldsymbol{q}}$. The
evaluation of these coefficients, as well as to obtain the eigenenergies,
requires a self-consistent calculation.

The analytical expressions obtained for the mean-field coefficients, Eqs.
(11)-(12) constitute the main result of our work. To illustrate the effects of
interacting thermal magnons on the AF magnon spectrum and the polarization of
magnons, we analyze both uniaxial (${\cal K}_{x}=0$) and biaxial (${\cal
K}_{x}\neq 0$) AF systems with hexagonal and square lattice structures that
are two common lattice geometries.

Figure 1: Temperature dependence of degenerate magnon eigenenergies of a
uniaxial AF system with (a) 2D square and (b) 2D hexagonal lattices as a
function of dimensionless wavevector $q_{x}/q_{m}$. Dashed curve lines
represent the bare magnon dispersion given by Eq. (4) with
$\mathcal{K}_{x}=0$. At the inset of panels (a) and (b), the AF resonance
modes (renormalized magnon bandgap) as a function of temperature are
presented. The mean-field coefficients (in units of meV) , ${\cal C}$ and
${\cal D}^{\prime}$, are shown for different temperatures as a function of
wavevector for square and hexagonal lattices at panels (c) and (d),
respectively. We use these typical parameters: $J=1$ [meV],
$\mathcal{K}_{z}=0.1$ [meV], $S=5/2$ and $a=1$ [nm].

##  IV Nonlinear spectral shift of magnons

In this section, we investigate the effect of nonlinear magnon interactions on
magnon dispersion, magnon bandgaps, and AF resonance modes of uniaxial and
biaxial antiferromagnetic insulators.

###  IV.1 Uniaxial AF systems

First, for the sake of completeness, we consider an AF system with easy-axis
anisotropy along the $z$-direction, modeled by the Hamiltonian in Eq. (1) when
$\mathcal{K}_{x}=0$. Nonlinear magnon interactions in uniaxial AF systems have
been previously treated in detail to unveil their role in the lifetime of
magnons Harris  _et al._ (1971); Bayrakci  _et al._ (2013), hydrodynamic
regime Harris  _et al._ (1971) and bulk spin transport Shen (2019); Troncoso
_et al._ (2020).

As we have already mentioned, in the uniaxial AF system and in the absence of
a magnetic field and DMIs, the two AF magnon eigenmodes are degenerate and
given by Eq. (4) with $\mathcal{K}_{x}=0$, and their corresponding eigenstates
are circularly polarized magnons. On the other hand, the effective interacting
Hamiltonian of magnons in the uniaxial AF case, in a basis that diagonalizes
the bare part and within the mean-field approximation, see Appendix A.1; is
given by,

$\displaystyle{H}^{\mathrm{u}(4)}_{\text{MF}}=\sum_{\boldsymbol{q}}\left[{\cal
C}_{\boldsymbol{q}}\alpha^{\mathrm{u}\dagger}_{\boldsymbol{q}}\alpha^{\mathrm{u}}_{\boldsymbol{q}}+{\cal
C}^{\prime}_{\boldsymbol{q}}\beta^{\mathrm{u}{\dagger}}_{\boldsymbol{q}}\beta^{\mathrm{u}}_{\boldsymbol{q}}+{\cal
D}^{\prime}_{\boldsymbol{q}}\alpha^{\mathrm{u}}_{\boldsymbol{q}}\beta^{\mathrm{u}}_{-\boldsymbol{q}}+\text{H.c.}\right],$
(13)

where the superscript $"\mathrm{u}"$ refers to the uniaxial case. The mean-
field coefficients ${\cal C}_{\boldsymbol{q}}$ and ${\cal
C}^{\prime}_{\boldsymbol{q}}$ renormalize two noninteracting magnon modes,
described by the linear spin-wave theory $H^{(2)}_{\text{m}}$, Eq. (3). In the
absence of magnetic fields and DMIs ${\cal C}_{\boldsymbol{q}}={\cal
C}^{\prime}_{\boldsymbol{q}}$.

The magnon nonconserving term with strength ${\cal
D}^{\prime}_{\boldsymbol{q}}$, that represents interband magnon scatterings
with opposite linear momenta, breaks the local $U(1)$ symmetry associated with
the magnon number conservation. This process is equivalent to the excitation
and annihilation of a magnon pair with zero net linear and spin angular
momenta. We later show how this term leads to the thermal squeezing of magnon
eigenmodes.

The temperature-dependent mean-field coefficients are given by,

$\displaystyle{\cal C}_{\boldsymbol{q}}$
$\displaystyle=\frac{1}{2}[A_{1}(|u^{\mathrm{u}}_{\boldsymbol{q}}|^{2}+|v^{\mathrm{u}}_{\boldsymbol{q}}|^{2})-(A_{2,\boldsymbol{q}}\,\bar{u}^{\mathrm{u}}_{\boldsymbol{q}}\bar{v}^{\mathrm{u}}_{\boldsymbol{q}}+\bar{A}_{2,\boldsymbol{q}}\,u^{\mathrm{u}}_{\boldsymbol{q}}v^{\mathrm{u}}_{\boldsymbol{q}})],$
(14a) $\displaystyle{\cal C}^{\prime}_{\boldsymbol{q}}$
$\displaystyle=\frac{1}{2}[A_{1}(|u^{\mathrm{u}}_{\boldsymbol{q}}|^{2}+|v^{\mathrm{u}}_{\boldsymbol{q}}|^{2})-(A_{2,-\boldsymbol{q}}\,\bar{u}^{\mathrm{u}}_{\boldsymbol{q}}\bar{v}^{\mathrm{u}}_{\boldsymbol{q}}+\bar{A}_{2,-\boldsymbol{q}}\,u^{\mathrm{u}}_{\boldsymbol{q}}v^{\mathrm{u}}_{\boldsymbol{q}})],$
(14b) $\displaystyle{\cal D}^{\prime}_{\boldsymbol{q}}$
$\displaystyle=-2A_{1}\bar{u}^{\mathrm{u}}_{\boldsymbol{q}}v^{\mathrm{u}}_{\boldsymbol{q}}+[\bar{A}_{2,\boldsymbol{q}}(v^{\mathrm{u}}_{\boldsymbol{q}})^{2}+A_{2,\boldsymbol{q}}(\bar{u}^{\mathrm{u}}_{\boldsymbol{q}})^{2}],$
(14c)

the Bogoliubov parameters of the uniaxial case,
$u^{\mathrm{u}}_{\boldsymbol{q}}$ and $v^{\mathrm{u}}_{\boldsymbol{q}}$, and
temperature-dependent mean-field coefficients $A_{1}$ and $A_{2}$ are
presented in Appendix B. Even at zero temperature, these coefficients are
finite. Therefore, we can generally say that there is always a finite
nonlinear quantum correction to the bare magnon dispersion in AF systems.

The resulting total magnon Hamiltonian of a uniaxial AF system in the mean-
field approximation is
$H^{\mathrm{u}}={H}^{\mathrm{u}(2)}_{\text{m}}+{H}^{\mathrm{u}(4)}_{\text{MF}}$.
Due to the presence of interband processes, ${\cal D}^{\prime}$, $\alpha$ and
$\beta$ modes are coupled and the Hamiltonian $H^{\mathrm{u}}$ is no longer
diagonal in the $\\{\alpha,\beta\\}$ basis. To diagonalize $H^{\mathrm{u}}$,
we should introduce a new diagonal basis of operators,
$\hat{\mu}_{\boldsymbol{q}}$ and $\hat{\nu}_{\boldsymbol{q}}$, and the
Hamiltonian in this new basis reads
$H^{\mathrm{u}}=\sum_{\boldsymbol{q}}\left(\varepsilon^{\mu}_{\boldsymbol{q}}\hat{\mu}^{\dagger}_{\boldsymbol{q}}\hat{\mu}_{\boldsymbol{q}}+\varepsilon^{\nu}_{\boldsymbol{q}}\hat{\nu}^{\dagger}_{\boldsymbol{q}}\hat{\nu}_{\boldsymbol{q}}\right)$.
In the absence of magnetic fields and DMIs, the effective eigenmodes of an
interacting uniaxial AF system remains degenerate,
$\varepsilon^{\mu}_{\boldsymbol{q}}=\varepsilon^{\nu}_{\boldsymbol{q}}\equiv\varepsilon_{\boldsymbol{q}}$,
where

$\displaystyle\varepsilon_{\boldsymbol{q}}=\sqrt{\left(\epsilon^{\text{u}}_{\boldsymbol{q}}+{\cal
C}_{\boldsymbol{q}}\right)^{2}-|{\cal D}^{\prime}_{{\boldsymbol{q}}}|^{2}},$
(15)

with
$\epsilon^{\text{u}}_{\boldsymbol{q}}=zJS\sqrt{\big{(}1+2\mathcal{K}_{z}/zJ\big{)}^{2}-|\gamma_{\boldsymbol{q}}|^{2}}$
denotes the dispersion relation of noninteracting magnons in the uniaxial
case, i.e., Eq. (4) with $\mathcal{K}_{x}=0$.

The dispersion in Eq. (15) is valid for two-sublattice AF systems with
arbitrary dimensions and lattice structures at finite temperature. Here, we
consider two distinct geometries of spin systems which are quite common, i.e.,
square and hexagonal lattices. The evaluation of mean-field coefficients, Eqs.
(14a)-(14c), and the eigenenergy $\varepsilon_{\boldsymbol{q}}$, is done self-
consistently and involves the calculation of $A_{1}$ and $A_{2}$ for various
temperatures.

In Fig. 1, we display the magnon eigenenergies and mean-field parameters as a
function of the wavevector for different temperatures and for two lattice
structures. As we have already discussed, the eigenenergies do not reduce to
the eigenenergy of noninteracting magnons
$\epsilon^{\text{u}}_{\boldsymbol{q}}$ at zero temperature (dashed line at
Fig. 1(a) and (b)), thanks to the existence of nonlinear quantum corrections.
At finite temperature, we observe that magnon interactions lower the energy of
long wavelength magnons. In particular, from the inset of panels (a) and (b),
we see that AF resonance modes are decreased by increasing the temperature.
This figure shows that nonlinear magnon interactions are more pronounced in AF
systems with hexagonal lattice than their square lattice counterparts.

The mean-field parameters ${\cal C}_{\boldsymbol{q}}$ and ${\cal
D}^{\prime}_{\boldsymbol{q}}$, plotted at Fig. 1(c) and (d) for each lattice,
are both peaked at the zero wavevector. The intraband coefficient ${\cal
C}_{\boldsymbol{q}}$ may change the sign by varying the temperature and
wavevectore while the interband coefficient ${\cal
D}^{\prime}_{\boldsymbol{q}}$ is always positive in our formalism.

Figure 2: Temperature dependence of magnon eigenenergies of a biaxial AF
system in a 2D square lattice. In panel (a), the eigenenergies for both $\mu$
and $\nu$ modes as a function of the wavevector are plotted. Dashed lines
represent the noninteracting eigenenergies,
$\epsilon^{\alpha}_{\boldsymbol{q}}$ and $\epsilon^{\beta}_{\boldsymbol{q}}$,
see Eq. (4). The AF resonance modes, $\epsilon^{\mu}_{0}$ and
$\epsilon^{\nu}_{0}$, are shown in their respective panels, as a function of
temperature and different hard-axis anisotropy strengths. In panel (b), the
mean-field coefficients (in meV units) ${\cal C}_{\alpha\alpha},{\cal
C}^{\prime}_{\alpha\alpha},{\cal C}_{\beta\beta}$ and ${\cal
C}^{\prime}_{\beta\beta}$, are plotted for various temperatures as a function
of wavevector. The parameters employed are $J=1$ [meV],
$\mathcal{K}_{x}=2\mathcal{K}_{z}=0.2$ [meV], $S=2.5$ and $a=1$ [nm].

###  IV.2 Biaxial AF systems

In this subsection, we focus on the AF systems with biaxial anisotropy.
Previously, we have shown that the degeneracy of two noninteracting magnon
modes is lifted in this case, see Eq. (4), and the corresponding eigenstates
are elliptically polarized magnons Rezende  _et al._ (2019). The total
interacting magnon Hamiltonian within the mean-field approximation is
$H=H^{(2)}_{\text{m}}+H^{(4)}_{\text{MF}}$. The mean-field Hamiltonian (10)
can be partitioned into the intra- and interband contribution as
$H^{(4)}_{\text{MF}}=H^{(4)}_{\text{intra}}+H^{(4)}_{\text{inter}}$, with

$\displaystyle{H}_{\text{intra}}^{(4)}=$
$\displaystyle\sum_{\boldsymbol{q}}\Bigl{[}{\cal
C}_{\alpha\alpha,{\boldsymbol{q}}}\alpha^{\dagger}_{\boldsymbol{q}}\alpha_{\boldsymbol{q}}+{\cal
C}_{\beta\beta,{\boldsymbol{q}}}\beta^{\dagger}_{\boldsymbol{q}}\beta_{\boldsymbol{q}}$
$\displaystyle+{\cal
C}^{\prime}_{\alpha\alpha,{\boldsymbol{q}}}\alpha_{\boldsymbol{q}}\alpha_{-{\boldsymbol{q}}}+{\cal
C}^{\prime}_{\beta\beta,{\boldsymbol{q}}}\beta_{\boldsymbol{q}}\beta_{-{\boldsymbol{q}}}+\text{H.c.}\Bigr{]},$
(16a) $\displaystyle{H}_{\text{inter}}^{(4)}=$
$\displaystyle\sum_{\boldsymbol{q}}\left[{\cal
D}_{\alpha\beta,{\boldsymbol{q}}}\alpha^{\dagger}_{\boldsymbol{q}}\beta_{\boldsymbol{q}}+{\cal
D}^{\prime}_{\alpha\beta,{\boldsymbol{q}}}\alpha_{\boldsymbol{q}}\beta_{-{\boldsymbol{q}}}+\text{H.c.}\right].$
(16b)

The intraband contribution of nonlinear magnon scattering consists of diagonal
elements, in the $\alpha$ and $\beta$ operators representation, with
coefficients ${\cal C}_{\chi\chi}$ and the off-diagonal elements with
coefficients ${\cal C}^{\prime}_{\chi\chi}$. ${\cal
C}_{\chi\chi=\alpha\alpha(\beta\beta)}$ renormalizes noninteracting magnon
dispersions, obtained within the linear spin-wave theory $H^{(2)}_{\text{m}}$,
Eq. (3), while ${\cal C}^{\prime}_{\chi\chi=\alpha\alpha(\beta\beta)}$
characterize magnon scatterings between states with opposite momenta, $\bm{q}$
and $-\bm{q}$, inside a particular magnon band. The later is equivalent to the
annihilation and excitation of an (a) $\alpha(\beta)-$magnon pair with zero
net linear momentum and spin angular momentum $+2\hbar(-2\hbar)$. The
interband nonlinear magnon scatterings are parametrized by off-diagonal terms,
${\cal D}_{\alpha\beta}$ and ${\cal D}^{\prime}_{\alpha\beta}$, corresponding
to scattering processes of magnons between $\alpha$ and $\beta$ bands. The
total interacting Hamiltonian, $H$, is no longer diagonalized in the basis of
$\alpha$ and $\beta$ operators. Thus, we do another canonical transformation
and find the diagonal basis in terms of new bosonic operators, say $\mu$ and
$\nu$. The expressions for mean-field parameters are presented in the Appendix
B.2, see Eqs. (27a)-(27d). In the biaxial case, similar to the uniaxial one,
the temperature-dependent mean-field coefficients are finite even at $T=0$.
This nonlinear quantum corrections to magnon dispersion scales as $1/N$, see
the Appendix B.2, and therefore are negligible for large systems. The
nonlinear quantum corrections, in addition to linear quantum spin
fluctuations, could play a relevant role in the quantum spin-transfer torques
and spin shot noise phenomena Qaiumzadeh and Brataas (2018). The importance of
this correction is, however, conditioned to small systems, which is a size
limit where our mean-field approach might breakdown.

Figure 3: Temperature dependence of magnon eigenenergies of a biaxial AF
system in a 2D hexagonal lattice. In panels (a) and (b), the energies
$\tilde{\epsilon}^{\mu}_{q}$ and $\tilde{\epsilon}^{\nu}_{q}$ are plotted as a
function of the wavevector. Dashed lines represent the respective
noninteracting energies, see Eq. (4). The mean-field coefficients (in meV
units) ${\cal C}_{\alpha\alpha},{\cal C}^{\prime}_{\alpha\alpha},{\cal
C}_{\beta\beta},{\cal C}^{\prime}_{\beta\beta},{\cal D}_{\alpha\beta}$ and
${\cal D}^{\prime}_{\alpha\beta}$, are shown as a function of the wavevector
in panels (c) and (d) with $k_{B}T=40[\text{meV}]$ and
$K_{x}=0.2[\text{meV}]$. Note that the real parts of coefficients ${\cal
D}_{\alpha\beta}$ and ${\cal D}^{\prime}_{\alpha\beta}$ are zero, and thus
only imaginary parts are displayed. The parameters employed are the same as in
Fig. 2.

We find that in square lattices ${\cal D}_{\alpha\beta}={\cal
D}^{\prime}_{\alpha\beta}=0$, hence, the interband scattering processes are
absent. The eigenenergies for the two magnon bands in a square lattice are
given by,

$\displaystyle\epsilon^{\mu}_{\boldsymbol{q}}$
$\displaystyle=\sqrt{(\epsilon^{\alpha}_{{\boldsymbol{q}}}+{\cal
C}_{\alpha\alpha,{\boldsymbol{q}}})^{2}-|{\cal
C}^{\prime}_{\alpha\alpha,{\boldsymbol{q}}}|^{2}},$ (17a)
$\displaystyle\epsilon^{\nu}_{\boldsymbol{q}}$
$\displaystyle=\sqrt{(\epsilon^{\beta}_{{\boldsymbol{q}}}+{\cal
C}_{\beta\beta,{\boldsymbol{q}}})^{2}-|{\cal
C}^{\prime}_{\beta\beta,{\boldsymbol{q}}}|^{2}}.$ (17b)

In Fig. 2(a), we show the magnon eigenenergy of a 2D square lattice as a
function of the wavevector for different temperatures. The situation is
similar in 3D AF systems, see the Appendix B.3. The noninteracting
eigenenergies, $\epsilon^{\alpha}_{\boldsymbol{q}}$ and
$\epsilon^{\beta}_{\boldsymbol{q}}$, in the upper- and lower-panel are
respectively represented by dashed lines. The effect of nonlinear interactions
is in particular evidenced by the temperature dependence of the AF resonance
modes, $\epsilon^{\mu}_{0}$ and $\epsilon^{\nu}_{0}$, shown in the respective
insets of panel (a). These figures show that the energy of one magnon band
$\epsilon^{\nu}_{0}$, decreases with temperature, while another band
$\epsilon^{\mu}_{0}$ has a weak temperature dependence for different strengths
of the hard-axis anisotropy.

Quite contrary, we find that in hexagonal lattices, the interband mean-field
coefficients ${\cal D}_{\alpha\beta}$ and ${\cal D}^{\prime}_{\alpha\beta}$
are nonzero. The eigenenergies for the magnonic $\mu$\- and $\nu$-eigenmodes,
denoted by $\tilde{\epsilon}^{\mu}_{\boldsymbol{q}}$ and
$\tilde{\epsilon}^{\nu}_{\boldsymbol{q}}$, respectively are shown in Fig. 3.
Their analytical expressions are displayed in Eq. (B.3) at Appendix B.3.

Generally, in biaxial cases, the degeneracy of two magnon bands are broken and
one band (the $\nu$ band) has lower energy than other one (the $\mu$ band).
Therefore, increasing temperature leads to an increase in the population of
the lowest magnon band, the $\nu$ band, and consequently these thermally
excited magnons renormaize the the $\nu$ band more strongly than the $\mu$
band. This can be readily seen in Figs. 2 and 3.

##  V Thermal Squeezed Magnons

Now, in this section, we study the effect of nonlinear thermal magnons on the
polarization of magnon eigenmodes. In general, magnon polarization is an
essential ingredient in magnonic-based information technology Qaiumzadeh  _et
al._ (2018); Jin  _et al._ (2021); Barman  _et al._ (2021). Therefore, it is
important to understand the effect of interactions on the polarization of
magnon eigenmodes. As we have mentioned earlier, two-sublattice AF systems
have two magnon eigenmodes. At long-wavelength limit, i.e., close to the
magnetic Brillouin zone center, the polarization of magnons is determined by
magnetic anisotropy, while at short-wavelength, i.e., close to the magnetic
Brillouin zone boundaries, the polarization of the modes are governed by the
Heisenberg exchange interactions. Therefore, in our model with isotropic
Heisenberg exchange interaction in the presence of either uniaxial or biaxial
magnetic anisotropies, high energy magnon eigenmodes close to the BZ magnetic
boundaries are circularly polarized with an ellipticity $e=\langle
S^{x}_{i}\rangle_{\rm{max}}/\langle S^{y}_{i}\rangle_{\rm{max}}=1$, see the
inset of Fig. 4 and right panel of Fig. 5. On the other hand, in the previous
section, we have shown that nonlinear thermal interactions mostly modify the
low-energy magnon excitations. Therefore, we expect interactions mainly affect
on the polarization of the anisotropy-dominated low-energy magnons.

It has been shown, both theoretically and experimentally, that in quantum AF
systems with negligible magnetic anisotropy and at low temperature, there are
magnon self-squeezing states because of AF quantum fluctuations and
interaction between AF sublattices Peng (2001); Zhao  _et al._ (2004, 2006);
Bossini  _et al._ (2019); Wuhrer  _et al._ (2022). Magnon squeezed states and
quantum spin entanglement are related, and thus they are interesting for
applications in quantum computing. In the presence of anisotropy and finite
temperature, magnon squeezing is reduced, but here we show that thermal
squeezing becomes dominant at finite temperature Cheng  _et al._ (2008). This
thermal squeezing in uniaxial AF systems, arising from an effective attractive
magnon interaction and breaking of the local U(1) symmetry, see Eq. (13) and
discussions below that.

Figure 4: Temperature-dependent ellipticity of AF resonant modes of a uniaxial
AF square lattice for different easy-axis anisotropy ($K_{z}$). The two AF
modes are degenerate with the same ellipticity $e_{\rm{uni}}$. The inset shows
ellepticity vs. dimensionless wavenumber for different temperatures,
$k_{B}T=\\{0.1,10,20,30\\}$ [meV].

The resulting temperature-dependent magnonic modes exhibit anisotropic spatial
oscillations. This property is evidenced by the ellipticity, which is
determinedRezende  _et al._ (2019) from the Bogolioubov transformation, (Eqs.
19 and 21), and diagonalization matrix of the Hamiltonian
$H={H}^{(2)}_{\text{m}}+{H}^{(4)}_{\text{MF}}$, that in turn relates the spin
components and diagonal basis, $\hat{\mu}$ and $\hat{\nu}$, of operators. The
result is obtained for the uniaxial ($e_{\rm{uni}}$) and biaxial
($e_{\rm{bi}}$) AF square-lattice systems, and presented at Figs. 4 and 5 as a
function of temperature for different magnetic anisotropy strengths. Similar
results can also be found in hexagonal lattices.

Figure 5: Ellipticity of two AF resonance modes as a function of temperature
(left panel) and dimensionless wavenumber (right panel) in a biaxial AF square
lattice. The temperature-dependent ellipticity of each AF resonance mode,
$e^{\mu}_{\rm{bi}}$ and $e^{\nu}_{\rm{bi}}$, are represented at left panels by
blue and red curves, respectively. The result is shown for various hard-axis
anisotropies, $K_{x}$, with a fixed easy axis anisotropy $K_{z}=0.1$. Right
panels show the wavenumber dependence of $e^{\mu}_{\rm{bi}}$ and
$e^{\nu}_{\rm{bi}}$ for different temperatures, $k_{B}T=\\{0,0.1,10,20,30\\}$
[meV] around the center of the first BZ. Black dashed line represents the
ellipticity of noninteracting magnons. The inset in the right panel shows the
behaviour of ellipticity in the whole BZ, where two magnon modes become
circularly polarized close to the first BZ boundaries.

In Fig. 4, we plot the ellipticity $e_{\rm{uni}}$ for the degenerate magnon
modes in uniaxial AF systems. The result is identical for each magnon mode,
which is shown for various easy-axis anisotropy values. Interestingly, the
eigenmode oscillations become squeezed along $y$-axis when temperature
increases. This effect is a direct consequence of U(1) symmetry breaking
arising from nonlinear thermal interactions, see Eq. (13), and thus is termed
thermal squeezed magnonsCheng  _et al._ (2008), unlike its quantum counterpart
defined at $T=0$ Peng (2001); Zhao  _et al._ (2004, 2006); Bossini  _et al._
(2019). In the biaxial AF system the ellipticity, $e_{\rm{bi}}$, of $\mu$\-
and $\nu$-magnonic eigenmode, are plotted at Fig. 5 as a function of
temperature and momentum. In this case, the eigenmodes are nondegenerate and
they are elliptical even at zero temperature along different directions. The
ellipticities, $e^{\mu}_{\rm{bi}}$ and $e^{\nu}_{\rm{bi}}$, at
${\boldsymbol{q}}={\boldsymbol{0}}$ are represented (left panels) by blue and
red curves, respectively. The magnonic $\mu$($\nu$)-mode becomes more (less)
elliptical, being $e^{\mu}_{\rm{bi}}>1$ and $e^{\nu}_{\rm{bi}}<1$, when
temperature increases but in an anisotropic way. Similarly, the dependence
with momentum is displayed at right panels. As temperature is increased, both
$e^{\mu}_{\rm{bi}}$ and $e^{\nu}_{\rm{bi}}$, slightly deviate from the
ellipticity for noninteracting magnons (i.e., at $T=0$ and depicted by a black
dashed line). Different rate in the change of ellipticity of the two modes vs.
temperature is related to the fact that the $\nu$ mode has lower energy and
this band becomes more populated than the $\mu$ mode; therefore the nonlinear
effects are stronger in this lower band, see Fig. 2.

##  VI Summary and Conclusion

We have developed a self-consistent Hartree-Fock mean-field formalism to
compute nonlinear magnon interactions in two-sublattice biaxial AF systems at
finite temperature. We have shown that thermal interactions, on the one hand,
reduce magnon bandgaps and thus AF resonance modes. Magnon bandgap plays an
important role in the stability of long-range magnetic order in 2D systems. On
the other hand, we have investigated the effect of nonlinear thermal magnon
interactions on the polarization of low-energy magnons. We have shown that in
the uniaxial AF case, magnon modes become squeezed because of the attractive
nature of nonlinear magnon interactions and the breaking of local U(1)
symmetry. In the present paper, we have used the truncated HP bosonization
transformation scheme to compute nonlinear magnons. The Dyson-Maleev
transformation Dyson (1956); Maleev (1958) also leads to the same effective
renormalized Hamiltonian within the Hartree-Fock approximation Liu (1992);
Mkhitaryan and Ke (2021). We have shown that long-wavelength magnons interact
more strongly than high-energy magnons. Even at zero temperature, there is a
finite nonlinear quantum correction to the bare magnon dispersion. The power
of this method is that it not only simplifies the interaction terms in the
magnon Hamiltonian but also does not impose a new constraint on the results.
The HP bosonization method is applied to find magnon interactions, and the
Hartree-Fock method, applied to find mean-field thermal magnon interactions,
is more accurate at low magnon numbers and large spin limits. Within this
formalism, it is easy to investigate the effect of nonlinear magnons in
various magnonic phenomena such as ultrafast magnon excitation and transport,
two-magnon dynamics, and magnon condensation.

###### Acknowledgements.

M. Sh. acknowledges Dr. Morteza Mohseni for fruitful discussions. M. Sh. and
J. F. acknowledge support from Carl Tryggers Stiftelse, and J. F. acknowledges
support from Vetenskapsrådet. This project has been supported by the Norwegian
Financial Mechanism Project No. 2019/34/H/ST3/00515, “2Dtronics”; and
partially by the Research Council of Norway through its Centres of Excellence
funding scheme, Project No. 262633, “QuSpin”.

## References

  * Marti  _et al._ (2014) X. Marti, I. Fina, C. Frontera, J. Liu, P. Wadley, Q. He, R. Paull, J. Clarkson, J. Kudrnovskỳ, I. Turek, _et al._ , Nature materials 13, 367 (2014).
  * Kriegner  _et al._ (2016) D. Kriegner, K. Vỳbornỳ, K. Olejník, H. Reichlová, V. Novák, X. Marti, J. Gazquez, V. Saidl, P. Němec, V. Volobuev, _et al._ , Nature communications 7, 1 (2016).
  * Lebrun  _et al._ (2018) R. Lebrun, A. Ross, S. A. Bender, A. Qaiumzadeh, L. Baldrati, J. Cramer, A. Brataas, R. A. Duine, and M. Kläui, Nature 561, 222–225 (2018).
  * Němec  _et al._ (2018) P. Němec, M. Fiebig, T. Kampfrath, and A. V. Kimel, Nature Physics 14, 229 (2018).
  * Geprags  _et al._ (2020) S. Geprags, M. Opel, J. Fischer, O. Gomonay, P. Schwenke, M. Althammer, H. Huebl, and R. Gross, Journal of Applied Physics 127, 243902 (2020).
  * Baltz  _et al._ (2018) V. Baltz, A. Manchon, M. Tsoi, T. Moriyama, T. Ono, and Y. Tserkovnyak, Rev. Mod. Phys. 90, 015005 (2018).
  * Qaiumzadeh  _et al._ (2018) A. Qaiumzadeh, L. A. Kristiansen, and A. Brataas, Phys. Rev. B 97, 020402 (2018).
  * Chumak  _et al._ (2014) A. V. Chumak, A. A. Serga, and B. Hillebrands, Nature communications 5, 1 (2014).
  * Bowlan  _et al._ (2018) P. Bowlan, S. A. Trugman, D. A. Yarotski, A. J. Taylor, and R. P. Prasankumar, Journal of Physics D: Applied Physics 51, 194003 (2018).
  * Jungwirth  _et al._ (2016) T. Jungwirth, X. Marti, P. Wadley, and J. Wunderlich, Nature nanotechnology 11, 231 (2016).
  * Proskurin  _et al._ (2017) I. Proskurin, R. L. Stamps, A. S. Ovchinnikov, and J.-i. Kishine, Phys. Rev. Lett. 119, 177202 (2017).
  * Chumak  _et al._ (2015) A. V. Chumak, V. I. Vasyuchka, A. A. Serga, and B. Hillebrands, Nature Physics 11, 453 (2015).
  * Gomonay and Loktev (2014) E. V. Gomonay and V. M. Loktev, Low Temperature Physics 40, 17 (2014).
  * Daniels  _et al._ (2015) M. W. Daniels, W. Guo, G. M. Stocks, D. Xiao, and J. Xiao, New Journal of Physics 17, 103039 (2015).
  * Khymyn  _et al._ (2016) R. Khymyn, I. Lisenkov, V. S. Tiberkevich, A. N. Slavin, and B. A. Ivanov, Phys. Rev. B 93, 224421 (2016).
  * Gomonay  _et al._ (2017) O. Gomonay, T. Jungwirth, and J. Sinova, physica status solidi (RRL) – Rapid Research Letters 11, 1700022 (2017).
  * Nishitani  _et al._ (2010) J. Nishitani, K. Kozuki, T. Nagashima, and M. Hangyo, Applied Physics Letters 96, 221906 (2010).
  * Nishitani  _et al._ (2012) J. Nishitani, T. Nagashima, and M. Hangyo, Phys. Rev. B 85, 174439 (2012).
  * Lebrun  _et al._ (2020) R. Lebrun, A. Ross, O. Gomonay, V. Baltz, U. Ebels, A.-L. Barra, A. Qaiumzadeh, A. Brataas, J. Sinova, and M. Kläui, Nature Communications 11 (2020), 10.1038/s41467-020-20155-7.
  * Rezende  _et al._ (2019) S. M. Rezende, A. Azevedo, and R. L. Rodríguez-Suárez, Journal of Applied Physics 126, 151101 (2019).
  * Bunkov and Volovik (2010) Y. M. Bunkov and G. E. Volovik, Journal of Physics: Condensed Matter 22, 164210 (2010).
  * Dzyapko  _et al._ (2017) O. Dzyapko, I. Lisenkov, P. Nowik-Boltyk, V. E. Demidov, S. O. Demokritov, B. Koene, A. Kirilyuk, T. Rasing, V. Tiberkevich, and A. Slavin, Phys. Rev. B 96, 064438 (2017).
  * Tupitsyn  _et al._ (2008) I. S. Tupitsyn, P. C. E. Stamp, and A. L. Burin, Phys. Rev. Lett. 100, 257202 (2008).
  * Rezende (2009) S. M. Rezende, Phys. Rev. B 79, 174411 (2009).
  * Troncoso  _et al._ (2015) R. E. Troncoso, C. Ulloa, F. Pesce, and A. S. Nunez, Phys. Rev. B 92, 224424 (2015).
  * Zakeri (2020) K. Zakeri, Journal of Physics: Condensed Matter 32, 363001 (2020).
  * Chumak  _et al._ (2017) A. V. Chumak, A. A. Serga, and B. Hillebrands, Journal of Physics D: Applied Physics 50, 244001 (2017).
  * Singh  _et al._ (2021) S. Singh, X. Wang, A. Kumar, A. Qaiumzadeh, P. Svedlindh, T. Tybell, and E. Wahlström, J. Magn. Magn. Mater. 539, 168376 (2021).
  * Kirilyuk  _et al._ (2010) A. Kirilyuk, A. V. Kimel, and T. Rasing, Rev. Mod. Phys. 82, 2731 (2010).
  * Wang  _et al._ (2020) K. Wang, X. Xu, Y. Cheng, M. Zhang, J.-S. Wang, H. Wang, and G. Zhang, Phys. Rev. B 102, 235434 (2020).
  * 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, Phys. Rev. Lett. 111, 017204 (2013).
  * Mohseni  _et al._ (2013) S. M. Mohseni, S. R. Sani, J. Persson, T. N. A. Nguyen, S. Chung, Y. Pogoryelov, P. K. Muduli, E. Iacocca, A. Eklund, R. K. Dumas, S. Bonetti, A. Deac, M. A. Hoefer, and J. Åkerman, Science 339, 1295 (2013).
  * Mohseni  _et al._ (2020) M. Mohseni, A. Qaiumzadeh, A. A. Serga, A. Brataas, B. Hillebrands, and P. Pirro, New Journal of Physics 22, 083080 (2020).
  * Arakawa (2019) N. Arakawa, Phys. Rev. B 99, 014405 (2019).
  * Giamarchi  _et al._ (2008) T. Giamarchi, C. Rüegg, and O. Tchernyshyov, Nature Physics 4, 198 (2008).
  * Mook  _et al._ (2021) A. Mook, K. Plekhanov, J. Klinovaja, and D. Loss, Phys. Rev. X 11, 021061 (2021).
  * Pershoguba  _et al._ (2018) S. S. Pershoguba, S. Banerjee, J. C. Lashley, J. Park, H. Ågren, G. Aeppli, and A. V. Balatsky, Phys. Rev. X 8, 011010 (2018).
  * Zhitomirsky and Chernyshev (2013) M. E. Zhitomirsky and A. L. Chernyshev, Rev. Mod. Phys. 85, 219 (2013).
  * Elliott and Thorpe (1969) R. Elliott and M. Thorpe, Journal of Physics C: Solid State Physics 2, 1630 (1969).
  * Liu (1992) B.-G. Liu, Journal of Physics: Condensed Matter 4, 8339 (1992).
  * Chernyshev and Zhitomirsky (2009) A. L. Chernyshev and M. E. Zhitomirsky, Phys. Rev. B 79, 144416 (2009).
  * McCreary  _et al._ (2020) A. McCreary, J. R. Simpson, T. T. Mai, R. D. McMichael, J. E. Douglas, N. Butch, C. Dennis, R. Valdés Aguilar, and A. R. Hight Walker, Phys. Rev. B 101, 064416 (2020).
  * Niyazi  _et al._ (2021) A. Niyazi, D. Geffroy, and J. Kunes, “Antiferromagnetic magnons and local anisotropy: dynamical mean-field study,” (2021), arXiv:2105.12468 [cond-mat.str-el] .
  * Mkhitaryan and Ke (2021) V. V. Mkhitaryan and L. Ke, Phys. Rev. B 104, 064435 (2021).
  * Li  _et al._ (2018) Z. Li, T. Cao, and S. G. Louie, Journal of Magnetism and Magnetic Materials 463, 28 (2018).
  * Mermin and Wagner (1966) N. D. Mermin and H. Wagner, Phys. Rev. Lett. 17, 1133 (1966).
  * Hohenberg (1967) P. C. Hohenberg, Phys. Rev. 158, 383 (1967).
  * Castilla and Chakravarty (1991) G. E. Castilla and S. Chakravarty, Phys. Rev. B 43, 13687 (1991).
  * Igarashi (1992) J.-i. Igarashi, Phys. Rev. B 46, 10763 (1992).
  * An  _et al._ (2001) J. An, Chang-DeGong, and H.-Q. Lin, J. Phys. Condens. Matter 13, 115 (2001).
  * Gitgeatpong  _et al._ (2017) G. Gitgeatpong, Y. Zhao, P. Piyawongwatthana, Y. Qiu, L. W. Harriger, N. P. Butch, T. J. Sato, and K. Matan, Phys. Rev. Lett. 119, 047201 (2017).
  * Qaiumzadeh  _et al._ (2017) A. Qaiumzadeh, H. Skarsvåg, C. Holmqvist, and A. Brataas, Phys. Rev. Lett. 118, 137201 (2017).
  * Holstein and Primakoff (1940) T. Holstein and H. Primakoff, Phys. Rev. 58, 1098 (1940).
  * König and Hucht (2021) J. König and A. Hucht, SciPost Phys. 10, 7 (2021).
  * Auerbach (1998) A. Auerbach,  _Interacting electrons and quantum magnetism_ (Springer-Verlag New York, 1998).
  * Harris  _et al._ (1971) A. B. Harris, D. Kumar, B. I. Halperin, and P. C. Hohenberg, Phys. Rev. B 3, 961 (1971).
  * Balucani  _et al._ (1980) U. Balucani, M. Pini, and V. Tognetti, Journal of Physics C: Solid State Physics 13, 2925 (1980).
  * Oguchi (1960) T. Oguchi, Phys. Rev. 117, 117 (1960).
  * Troncoso  _et al._ (2020) R. E. Troncoso, S. A. Bender, A. Brataas, and R. A. Duine, Phys. Rev. B 101, 054404 (2020).
  * Flebus (2019) B. Flebus, Phys. Rev. B 100, 064410 (2019).
  * Shen (2019) K. Shen, Phys. Rev. B 100, 094423 (2019).
  * Qaiumzadeh and Brataas (2018) A. Qaiumzadeh and A. Brataas, Phys. Rev. B 98, 220408 (2018).
  * Jin  _et al._ (2021) Z. Jin, C. Y. Meng, T. T. Liu, D. Y. Chen, Z. Fan, M. Zeng, X. B. Lu, X. S. Gao, M. H. Qin, and J.-M. Liu, Phys. Rev. B 104, 054419 (2021).
  * Barman  _et al._ (2021) A. Barman, G. Gubbiotti, S. Ladak, A. O. Adeyeye, M. Krawczyk, J. Gräfe, C. Adelmann, S. Cotofana, A. Naeemi, V. I. Vasyuchka, B. Hillebrands, S. A. Nikitov, H. Yu, D. Grundler, A. V. Sadovnikov, A. A. Grachev, S. E. Sheshukova, J.-Y. Duquesne, M. Marangolo, G. Csaba, W. Porod, V. E. Demidov, S. Urazhdin, S. O. Demokritov, E. Albisetti, D. Petti, R. Bertacco, H. Schultheiss, V. V. Kruglyak, V. D. Poimanov, S. Sahoo, J. Sinha, H. Yang, M. Münzenberg, T. Moriyama, S. Mizukami, P. Landeros, R. A. Gallardo, G. Carlotti, J.-V. Kim, R. L. Stamps, R. E. Camley, B. Rana, Y. Otani, W. Yu, T. Yu, G. E. W. Bauer, C. Back, G. S. Uhrig, O. V. Dobrovolskiy, B. Budinska, H. Qin, S. van Dijken, A. V. Chumak, A. Khitun, D. E. Nikonov, I. A. Young, B. W. Zingsem, and M. Winklhofer, J. Phys. Condens. Matter 33, 413001 (2021).
  * Peng (2001) F. Peng, Europhysics Letters (EPL) 54, 688 (2001).
  * Zhao  _et al._ (2004) J. Zhao, A. V. Bragas, D. J. Lockwood, and R. Merlin, Phys. Rev. Lett. 93, 107203 (2004).
  * Zhao  _et al._ (2006) J. Zhao, A. V. Bragas, R. Merlin, and D. J. Lockwood, Phys. Rev. B 73, 184434 (2006).
  * Bossini  _et al._ (2019) D. Bossini, S. Dal Conte, G. Cerullo, O. Gomonay, R. V. Pisarev, M. Borovsak, D. Mihailovic, J. Sinova, J. H. Mentink, T. Rasing, and A. V. Kimel, Phys. Rev. B 100, 024428 (2019).
  * Wuhrer  _et al._ (2022) D. Wuhrer, N. Rohling, and W. Belzig, Phys. Rev. B 105, 054406 (2022).
  * Cheng  _et al._ (2008) Z. Z. Cheng, B. Xu, Z. M. Li, and Z. Cheng, Eur. Phys. J. B 66, 289 (2008).
  * Dyson (1956) F. J. Dyson, Phys. Rev. 102, 1217 (1956).
  * Maleev (1958) S. V. Maleev, Soviet Physics JETP 6, 776 (1958).
  * Castro  _et al._ (2006) E. V. Castro, N. M. R. Peres, K. S. D. Beach, and A. W. Sandvik, Phys. Rev. B 73, 054422 (2006).
  * White  _et al._ (1965) R. M. White, M. Sparks, and I. Ortenburger, Phys. Rev. 139, A450 (1965).

##  Appendix A Bogoliubov Coefficients

In this Appendix, we present the elements of Bogoliubov transformation matrix
that diagonalizes the noninteracting bosonic Hamiltonian. The quadratic
noninteracting bosinic Hamiltonian, derived from the spin Hamiltonian (1), in
the momentum space is given by,

$\displaystyle\begin{split}H^{(2)}_{b}=&S\sum_{\boldsymbol{q}}\bigg{[}\big{(}zJ+\mathcal{K}_{x}+2\mathcal{K}_{z}\big{)}\big{(}a^{\dagger}_{\boldsymbol{q}}a_{\boldsymbol{q}}+b^{\dagger}_{\boldsymbol{q}}b_{\boldsymbol{q}}\big{)}+zJ\gamma_{-{\boldsymbol{q}}}a_{\boldsymbol{q}}b_{-{\boldsymbol{q}}}\\\
&+zJ\gamma_{{\boldsymbol{q}}}b^{\dagger}_{-{\boldsymbol{q}}}a^{\dagger}_{\boldsymbol{q}}+\frac{\mathcal{K}_{x}}{2}\big{(}a_{\boldsymbol{q}}a_{-{\boldsymbol{q}}}+b_{\boldsymbol{q}}b_{-{\boldsymbol{q}}}+a^{\dagger}_{\boldsymbol{q}}a^{\dagger}_{-{\boldsymbol{q}}}+b^{\dagger}_{\boldsymbol{q}}b^{\dagger}_{-{\boldsymbol{q}}}\big{)}\bigg{]}\end{split}$
(18)

where the following Fourier transform of operators is implemented
$a_{i}=\sqrt{{2}/{N}}\sum_{\boldsymbol{q}}e^{i{\boldsymbol{q}}\cdot\mathbf{r}_{i}}a_{\boldsymbol{q}}$
and
$b_{j}=\sqrt{{2}/{N}}\sum_{\boldsymbol{q}}e^{i{\boldsymbol{q}}\cdot\mathbf{r}_{j}}b_{\boldsymbol{q}}$,
with $\\{i,j\\}\in\\{\mathcal{A},\mathcal{B}\\}$. In order to explore four-
magnon scattering processes, one needs to find the magnon band structure. The
diagonalization of the Hamiltonian Eq. (18) is obtained through the canonical
Bogolioubov transformations Castro  _et al._ (2006) that maps the bosonic
operator $\\{a_{\boldsymbol{q}},b_{\boldsymbol{q}}\\}$ into the magnon
eigenbasis $\\{\alpha_{\boldsymbol{q}},\beta_{\boldsymbol{q}}\\}$. In the
following, we present the results for uniaxial and biaxial AF cases,
separately.

###  A.1 Uniaxial AF system

To diagonalize the noninteracting bosonic Hamiltonian (18) in the uniaxial
easy-axis limit ($\mathcal{K}_{x}=0$), we perform a canonical transformation
using the following $2\times 2$ Bogolioubov transformation, Shen (2019);
Rezende  _et al._ (2019),

$\displaystyle\begin{pmatrix}a_{\boldsymbol{q}}\\\
b^{\dagger}_{-{\boldsymbol{q}}}\end{pmatrix}=\begin{pmatrix}\bar{u}^{\mathrm{u}}_{\boldsymbol{q}}&-\bar{v}^{\mathrm{u}}_{\boldsymbol{q}}\\\
-v^{\mathrm{u}}_{\boldsymbol{q}}&u^{\mathrm{u}}_{\boldsymbol{q}}\end{pmatrix}\begin{pmatrix}\alpha^{\mathrm{u}}_{\boldsymbol{q}}\\\
\beta^{\mathrm{u}{\dagger}}_{\boldsymbol{q}}\end{pmatrix},$ (19)

where
$|u^{\mathrm{u}}_{\boldsymbol{q}}|^{2}-|v^{\mathrm{u}}_{\boldsymbol{q}}|^{2}=1$.
We define $u^{\mathrm{u}}_{\boldsymbol{q}}=\cosh(\theta_{\boldsymbol{q}}/2)$
and $v^{\mathrm{u}}_{\boldsymbol{q}}=\sinh(\theta_{\boldsymbol{q}}/2)$.
Inserting Eq. (19) into the Hamiltonian (18), we obtain the following
Bogoliubov coefficients for the uniaxial system,

$\displaystyle u^{\mathrm{u}}_{\boldsymbol{q}}$
$\displaystyle=\sqrt{\frac{zJS+2\mathcal{K}_{z}S+\epsilon^{u}_{\alpha,{\boldsymbol{q}}}}{2\epsilon^{u}_{\alpha,{\boldsymbol{q}}}}},$
(20a) $\displaystyle v^{\mathrm{u}}_{\boldsymbol{q}}$
$\displaystyle=\sqrt{\frac{zJS+2\mathcal{K}_{z}S-\epsilon^{u}_{\beta,{\boldsymbol{q}}}}{2\epsilon^{u}_{\beta,{\boldsymbol{q}}}}},$
(20b)

where the degenerate eigenenergies are
$\epsilon^{\mathrm{u}}_{\chi,{\boldsymbol{q}}}=S\sqrt{\big{(}zJ+2\mathcal{K}_{z}\big{)}^{2}-\big{(}zJ|\gamma_{\boldsymbol{q}}|\big{)}^{2}}$,
with $\chi\in\\{\alpha,\beta\\}$.

###  A.2 Biaxial AF system

In the case of biaxial AF systems, we use the following $4\times 4$
Bogolioubov transformations White  _et al._ (1965); Rezende  _et al._ (2019),

$\begin{split}\begin{pmatrix}a_{\boldsymbol{q}}\\\ b_{\boldsymbol{q}}\\\
a^{\dagger}_{-{\boldsymbol{q}}}\\\
b^{\dagger}_{-{\boldsymbol{q}}}\end{pmatrix}=\begin{pmatrix}u_{\boldsymbol{q}}&&v_{\boldsymbol{q}}&&x_{\boldsymbol{q}}&&w_{\boldsymbol{q}}\\\
-u_{\boldsymbol{q}}&&v_{\boldsymbol{q}}&&-x_{\boldsymbol{q}}&&w_{\boldsymbol{q}}\\\
\bar{x}_{\boldsymbol{q}}&&\bar{w}_{\boldsymbol{q}}&&\bar{u}_{\boldsymbol{q}}&&\bar{v}_{\boldsymbol{q}}\\\
-\bar{x}_{\boldsymbol{q}}&&\bar{w}_{\boldsymbol{q}}&&-\bar{u}_{\boldsymbol{q}}&&\bar{v}_{\boldsymbol{q}}\end{pmatrix}\begin{pmatrix}\alpha_{\boldsymbol{q}}\\\
\beta_{\boldsymbol{q}}\\\ \alpha^{\dagger}_{-{\boldsymbol{q}}}\\\
\beta^{\dagger}_{-{\boldsymbol{q}}}\end{pmatrix}\end{split}$ (21)

where the Bogoliubov coefficients are given by,

$\displaystyle u_{\boldsymbol{q}}$
$\displaystyle=\frac{1}{2}\sqrt{\frac{zJS+\mathcal{K}_{x}S+2\mathcal{K}_{z}S+\epsilon_{\alpha,{\boldsymbol{q}}}}{\epsilon_{\alpha,{\boldsymbol{q}}}}},$
(22a) $\displaystyle v_{\boldsymbol{q}}$
$\displaystyle=-\frac{1}{2}\sqrt{\frac{zJS+\mathcal{K}_{x}S+2\mathcal{K}_{z}S+\epsilon_{\beta,{\boldsymbol{q}}}}{\epsilon_{\beta,{\boldsymbol{q}}}}},$
(22b) $\displaystyle x_{\boldsymbol{q}}$
$\displaystyle=\frac{1}{2}\sqrt{\frac{zJS+\mathcal{K}_{x}S+2\mathcal{K}_{z}S-\epsilon_{\alpha,{\boldsymbol{q}}}}{\epsilon_{\alpha,{\boldsymbol{q}}}}},$
(22c) $\displaystyle w_{\boldsymbol{q}}$
$\displaystyle=\frac{1}{2}\sqrt{\frac{zJS+\mathcal{K}_{x}S+2\mathcal{K}_{z}S-\epsilon_{\beta,{\boldsymbol{q}}}}{\epsilon_{\beta,{\boldsymbol{q}}}}},$
(22d)

with
$\epsilon_{\chi,{\boldsymbol{q}}}=S\sqrt{\big{(}zJ+\mathcal{K}_{x}+2\mathcal{K}_{z}\big{)}^{2}-\big{(}zJ|\gamma_{\boldsymbol{q}}|\mp\mathcal{K}_{x}\big{)}^{2}}$
for magnon mode $\chi\in\\{\alpha,\beta\\}$. The quantum correction to the
ground-state energy is given by,

$\displaystyle\begin{split}E_{0}^{q}=S\sum_{\bf
q}\Big{\\{}2(zJ+\mathcal{K}_{x}+2\mathcal{K}_{z})(|x_{\bf q}|^{2}+|w_{\bf
q}|^{2})+2\,\text{Re}\big{[}zJ\gamma_{-\bf q}(-u_{\bf q}x_{\bf q}+v_{\bf
q}w_{\bf q})+\mathcal{K}_{x}(u_{\bf q}x_{\bf q}+v_{\bf q}w_{\bf
q})\big{]}\Big{\\}}.\end{split}$ (23)

##  Appendix B Mean-field Coefficients

In this Appendix, we present the treatment of the magnon-magnon Hamiltonian
$H^{(4)}_{b}$ within the mean-field approximation. We define the following
mean-field parametersBalucani  _et al._ (1980),

$\displaystyle m=\frac{2}{NS}\sum_{\boldsymbol{q}}\langle
a^{\dagger}_{\boldsymbol{q}}a_{\boldsymbol{q}}\rangle=\frac{2}{NS}\sum_{\boldsymbol{q}}\langle
b^{\dagger}_{\boldsymbol{q}}b_{\boldsymbol{q}}\rangle,$ (24a) $\displaystyle
m^{\prime}=\frac{2}{NS}\sum_{\boldsymbol{q}}\gamma_{-{\boldsymbol{q}}}\langle
a_{\boldsymbol{q}}b^{\dagger}_{\boldsymbol{q}}\rangle,$ (24b) $\displaystyle
p=\frac{2}{NS}\sum_{\boldsymbol{q}}\langle
a_{\boldsymbol{q}}a_{-{\boldsymbol{q}}}\rangle=\frac{2}{NS}\sum_{\boldsymbol{q}}\langle
b_{\boldsymbol{q}}b_{-{\boldsymbol{q}}}\rangle,$ (24c) $\displaystyle
p^{\prime}=\frac{2}{NS}\sum_{\boldsymbol{q}}\gamma_{-{\boldsymbol{q}}}\langle
a_{\boldsymbol{q}}b_{-{\boldsymbol{q}}}\rangle,$ (24d)

where $\langle...\rangle$ is for the thermal average. The mean-field parameter
$m$ denotes the number of bosonic excitations on each sublattice $\cal{A}$ and
$\cal{B}$, while $m^{\prime}$, $p$ and $p^{\prime}$ denote the interaction
between them. In the following, we evaluate these mean-field parameters for
both uniaxial and biaxial cases, separately.

###  B.1 Uniaxial AF system

In the uniaxial AF case, where the magnon basis reads as Eq. (19), the mean-
field parameters become,

$\displaystyle m$
$\displaystyle=\frac{2}{NS}\sum_{\boldsymbol{q}}\big{(}|u^{u}_{\boldsymbol{q}}|^{2}+|v^{u}_{\boldsymbol{q}}|^{2}\big{)}\,n_{\boldsymbol{q}}+|v^{u}_{\boldsymbol{q}}|^{2},$
(25a) $\displaystyle m^{\prime}$ $\displaystyle=p=0,$ (25b) $\displaystyle
p^{\prime}$
$\displaystyle=-\frac{2}{NS}\sum_{\boldsymbol{q}}\gamma_{-{\boldsymbol{q}}}u^{u}_{\boldsymbol{q}}v^{u}_{\boldsymbol{q}}\,(2n_{\boldsymbol{q}}+1).$
(25c)

The mean-field coefficients, introduced in the Hamiltonian (13) obey,

$\displaystyle{\cal C}_{\boldsymbol{q}}$
$\displaystyle=\frac{1}{2}[A_{1}(|u^{u}_{\boldsymbol{q}}|^{2}+|v^{u}_{\boldsymbol{q}}|^{2})-(A_{2,\boldsymbol{q}}\,\bar{u}^{u}_{\boldsymbol{q}}\bar{v}^{u}_{\boldsymbol{q}}+\bar{A}_{2,\boldsymbol{q}}\,u^{u}_{\boldsymbol{q}}v^{u}_{\boldsymbol{q}})],$
(26a) $\displaystyle{\cal C}^{\prime}_{\boldsymbol{q}}$
$\displaystyle=\frac{1}{2}[A_{1}(|u^{u}_{\boldsymbol{q}}|^{2}+|v^{u}_{\boldsymbol{q}}|^{2})-(A_{2,-\boldsymbol{q}}\,\bar{u}^{u}_{\boldsymbol{q}}\bar{v}^{u}_{\boldsymbol{q}}+\bar{A}_{2,-\boldsymbol{q}}\,u^{u}_{\boldsymbol{q}}v^{u}_{\boldsymbol{q}})],$
(26b) $\displaystyle{\cal D}^{\prime}_{\boldsymbol{q}}$
$\displaystyle=-2A_{1}\bar{u}^{u}_{\boldsymbol{q}}v^{u}_{\boldsymbol{q}}+(\bar{A}_{2,\boldsymbol{q}}(v^{u}_{\boldsymbol{q}})^{2}+A_{2,\boldsymbol{q}}(\bar{u}^{u}_{\boldsymbol{q}})^{2}).$
(26c)

where $A_{1}=-\big{(}zJ+4\mathcal{K}_{z}\big{)}Sm-
zJS(p^{\prime}+\bar{p}^{\prime})/2$ and
$A_{2}=-zJS(\bar{p}^{\prime}+m)\,\gamma_{-\boldsymbol{q}}$, and
$n_{\boldsymbol{q}}=(n_{\alpha,\boldsymbol{q}}+n_{\beta,\boldsymbol{q}})/2$.
In the uniaxial case and in the absence of magnetic fields and DMIs
$n_{\alpha,\boldsymbol{q}}=n_{\beta,\boldsymbol{q}}$.

###  B.2 Biaxial AF system

In this case ${\cal K}_{x}\neq 0$, where the magnon basis reads as Eq. (21),
the mean-field parameters become,

$\displaystyle m$
$\displaystyle=\frac{2}{NS}\sum_{\boldsymbol{q}}\left[\big{(}|u_{\boldsymbol{q}}|^{2}+|x_{\boldsymbol{q}}|^{2}\big{)}\,n_{\alpha,{\boldsymbol{q}}}+\big{(}|v_{\boldsymbol{q}}|^{2}+|w_{\boldsymbol{q}}|^{2}\big{)}\,n_{\beta,{\boldsymbol{q}}}+|x_{\boldsymbol{q}}|^{2}+|w_{\boldsymbol{q}}|^{2}\right],$
(27a) $\displaystyle p$
$\displaystyle=\frac{2}{NS}\sum_{\boldsymbol{q}}\left[u_{\boldsymbol{q}}x_{\boldsymbol{q}}\,\big{(}2n_{\alpha,{\boldsymbol{q}}}+1\big{)}+v_{\boldsymbol{q}}w_{\boldsymbol{q}}\,\big{(}2n_{\beta,{\boldsymbol{q}}}+1\big{)}\right],$
(27b) $\displaystyle m^{\prime}$
$\displaystyle=\frac{2}{NS}\sum_{\boldsymbol{q}}\gamma_{-{\boldsymbol{q}}}\left[-\big{(}|u_{\boldsymbol{q}}|^{2}+|x_{\boldsymbol{q}}|^{2}\big{)}\,n_{\alpha,{\boldsymbol{q}}}+\big{(}|v_{\boldsymbol{q}}|^{2}+|w_{\boldsymbol{q}}|^{2}\big{)}\,n_{\beta,{\boldsymbol{q}}}-|u_{\boldsymbol{q}}|^{2}+|v_{\boldsymbol{q}}|^{2}\right],$
(27c) $\displaystyle p^{\prime}$
$\displaystyle=\frac{2}{NS}\sum_{\boldsymbol{q}}\gamma_{-{\boldsymbol{q}}}\left[-u_{\boldsymbol{q}}x_{\boldsymbol{q}}\,\big{(}2n_{\alpha,{\boldsymbol{q}}}+1\big{)}+v_{\boldsymbol{q}}w_{\boldsymbol{q}}\,\big{(}2n_{\beta,{\boldsymbol{q}}}+1\big{)}\right],$
(27d)

which depend on the Bose-Einstein distribution function of magnons at each
mode, $n_{\alpha,{\boldsymbol{q}}}$ and $n_{\beta,{\boldsymbol{q}}}$, and the
Bogoliubov coefficients. In particular, note that the mean-field parameters
are finite even at zero temperature. In this regime, the magnon distributions
$n_{\chi,{\boldsymbol{q}}}\rightarrow 0$, and thus
$m=2\sum_{\boldsymbol{q}}(|x_{\boldsymbol{q}}|^{2}+|w_{\boldsymbol{q}}|^{2})/NS$,
$p=2\sum_{\boldsymbol{q}}(u_{\boldsymbol{q}}x_{\boldsymbol{q}}+v_{\boldsymbol{q}}w_{\boldsymbol{q}})/NS$,
$m^{\prime}=2\sum_{\boldsymbol{q}}\gamma_{-{\boldsymbol{q}}}(-|u_{\boldsymbol{q}}|^{2}+|v_{\boldsymbol{q}}|^{2})/NS$
and
$p^{\prime}=2\sum_{\boldsymbol{q}}\gamma_{-{\boldsymbol{q}}}(-u_{\boldsymbol{q}}x_{\boldsymbol{q}}+v_{\boldsymbol{q}}w_{\boldsymbol{q}})/NS$.
This results emphasize that magnon-magnon scattering introduces finite
corrections in the magnon dispersion even at zero temperature.
The mean-field coefficients, introduced in the Hamiltonian (10) obey,

$\displaystyle\cal{C}_{\alpha\alpha,{\boldsymbol{q}}}=$
$\displaystyle\frac{1}{2}\big{[}2A_{a}(|u_{\boldsymbol{q}}|^{2}+|x_{\boldsymbol{q}}|^{2})-(A^{\prime}_{ab,{\boldsymbol{q}}}+\bar{A}^{\prime}_{ab,{\boldsymbol{q}}})\,|u_{\boldsymbol{q}}|^{2}-(A^{\prime}_{ab,{-\boldsymbol{q}}}+\bar{A}^{\prime}_{ab,{-\boldsymbol{q}}})\,|x_{\boldsymbol{q}}|^{2}$
(28a)
$\displaystyle+(2A^{\prime}_{a}+2A^{\prime}_{b}-A_{ab,\boldsymbol{q}}-A_{ab,-\boldsymbol{q}})\,u_{\boldsymbol{q}}x_{\boldsymbol{q}}+(2\bar{A}^{\prime}_{a}+2\bar{A}^{\prime}_{b}-\bar{A}_{ab,\boldsymbol{q}}-\bar{A}_{ab,-\boldsymbol{q}})\,\bar{u}_{\boldsymbol{q}}\bar{x}_{\boldsymbol{q}}\big{]},$
$\displaystyle{\cal C}_{\beta\beta,{\boldsymbol{q}}}=$
$\displaystyle\frac{1}{2}\big{[}2A_{a}(|v_{\boldsymbol{q}}|^{2}+|w_{\boldsymbol{q}}|^{2})+(A^{\prime}_{ab,{\boldsymbol{q}}}+\bar{A}^{\prime}_{ab,{\boldsymbol{q}}})\,|v_{\boldsymbol{q}}|^{2}+(A^{\prime}_{ab,{-\boldsymbol{q}}}+\bar{A}^{\prime}_{ab,{-\boldsymbol{q}}})\,|w_{\boldsymbol{q}}|^{2}$
(28b)
$\displaystyle+(2A^{\prime}_{a}+2A^{\prime}_{b}+A_{ab,\boldsymbol{q}}+A_{ab,-\boldsymbol{q}})\,v_{\boldsymbol{q}}w_{\boldsymbol{q}}+(2\bar{A}^{\prime}_{a}+2\bar{A}^{\prime}_{b}+\bar{A}_{ab,\boldsymbol{q}}+\bar{A}_{ab,-\boldsymbol{q}})\,\bar{v}_{\boldsymbol{q}}\bar{w}_{\boldsymbol{q}}\big{]},$
$\displaystyle\cal{C}^{\prime}_{\alpha\alpha,{\boldsymbol{q}}}=$
$\displaystyle(2A_{a}-A^{\prime}_{ab,{\boldsymbol{q}}}-{\bar{A}^{\prime}}_{ab,{\boldsymbol{q}}})u_{\boldsymbol{q}}\bar{x}_{\boldsymbol{q}}+(A^{\prime}_{a}+A^{\prime}_{b}-A_{ab,{\boldsymbol{q}}})u^{2}_{\boldsymbol{q}}+({\bar{A}^{\prime}_{a}}+{\bar{A}^{\prime}_{b}}-\bar{A}_{ab,{\boldsymbol{q}}}){\bar{x}}^{2}_{\boldsymbol{q}},$
(28c) $\displaystyle\cal{C}^{\prime}_{\beta\beta,{\boldsymbol{q}}}=$
$\displaystyle(2A_{a}+A^{\prime}_{ab,{\boldsymbol{q}}}+{\bar{A}^{\prime}}_{ab,{\boldsymbol{q}}})v_{\boldsymbol{q}}\bar{w}_{\boldsymbol{q}}+(A^{\prime}_{a}+A^{\prime}_{b}+A_{ab,{\boldsymbol{q}}})v^{2}_{\boldsymbol{q}}+({\bar{A}^{\prime}_{a}}+{\bar{A}^{\prime}_{b}}+\bar{A}_{ab,{\boldsymbol{q}}}){\bar{w}}^{2}_{\boldsymbol{q}},$
(28d) $\displaystyle\cal{D}_{\alpha\beta,{\boldsymbol{q}}}=$
$\displaystyle(A_{ab,-{\boldsymbol{q}}}-A_{ab,{\boldsymbol{q}}})x_{\boldsymbol{q}}v_{\boldsymbol{q}}-(\bar{A}_{ab,-{\boldsymbol{q}}}-\bar{A}_{ab,{\boldsymbol{q}}})\bar{u}_{\boldsymbol{q}}\bar{w}_{\boldsymbol{q}}+\big{(}A^{\prime}_{ab,-{\boldsymbol{q}}}-{\bar{A}^{\prime}}_{ab,-{\boldsymbol{q}}}\big{)}x_{\boldsymbol{q}}\bar{w}_{\boldsymbol{q}}-\big{(}A^{\prime}_{ab,{\boldsymbol{q}}}-{\bar{A}^{\prime}}_{ab,{\boldsymbol{q}}}\big{)}v_{\boldsymbol{q}}\bar{u}_{\boldsymbol{q}}$
(28e)
$\displaystyle+2(A^{\prime}_{a}-A^{\prime}_{b})v_{\boldsymbol{q}}x_{\boldsymbol{q}}+2(\bar{A}^{\prime}_{a}-\bar{A}^{\prime}_{b})\bar{u}_{\boldsymbol{q}}\bar{w}_{\boldsymbol{q}},$
$\displaystyle\cal{D}^{\prime}_{\alpha\beta,{\boldsymbol{q}}}=$
$\displaystyle(\bar{A}_{ab,-{\boldsymbol{q}}}-\bar{A}_{ab,{\boldsymbol{q}}})\bar{x}_{\boldsymbol{q}}\bar{w}_{\boldsymbol{q}}+(A_{ab,{\boldsymbol{q}}}-A_{ab,-{\boldsymbol{q}}})u_{\boldsymbol{q}}v_{\boldsymbol{q}}+\big{(}{\bar{A}^{\prime}_{ab,-{\boldsymbol{q}}}}-A^{\prime}_{ab,-{\boldsymbol{q}}}\big{)}\bar{x}_{\boldsymbol{q}}v_{\boldsymbol{q}}+\big{(}A^{\prime}_{ab,{\boldsymbol{q}}}-{\bar{A}^{\prime}}_{ab,{\boldsymbol{q}}}\big{)}u_{\boldsymbol{q}}\bar{w}_{\boldsymbol{q}}$
(28f)
$\displaystyle+2(A^{\prime}_{a}-A^{\prime}_{b})u_{\boldsymbol{q}}v_{\boldsymbol{q}}+2(\bar{A}^{\prime}_{a}-\bar{A}^{\prime}_{b})\bar{x}_{\boldsymbol{q}}\bar{w}_{\boldsymbol{q}},$

where the following definitions have been introduced

$\displaystyle A_{a}$
$\displaystyle=-S\left[zJ\left(m+\frac{p^{\prime}+\bar{p}^{\prime}}{2}\right)+4m\left(\mathcal{K}_{z}+\frac{\mathcal{K}_{x}}{2}\right)-\frac{3}{4}\left(p+\bar{p}\right)\mathcal{K}_{x}\right]$
(29a) $\displaystyle A^{\prime}_{a}$
$\displaystyle=-S\left[zJ\bar{m}^{\prime}+\bar{p}\left(\mathcal{K}_{z}+\frac{\mathcal{K}_{x}}{2}\right)-\frac{3}{4}m\mathcal{K}_{x}\right],$
(29b) $\displaystyle A^{\prime}_{b}$
$\displaystyle=-S\left[zJ{m^{\prime}}+\bar{p}\left(\mathcal{K}_{z}+\frac{\mathcal{K}_{x}}{2}\right)-\frac{3}{4}m\mathcal{K}_{x}\right],$
(29c) $\displaystyle A_{ab,{\boldsymbol{q}}}$
$\displaystyle=-zJS\,\gamma_{-{\boldsymbol{q}}}\,\left(\bar{p}^{\prime}+m\right)$
(29d) $\displaystyle A^{\prime}_{ab,{\boldsymbol{q}}}$
$\displaystyle=-zJS\,\gamma_{-{\boldsymbol{q}}}\,\left(\bar{m}^{\prime}+\frac{p+\bar{p}}{4}\right).$
(29e)

In the cases evaluated in the main text, i.e., in the absence of magnetic
field and DMIs, the Bogoliubov coefficients
$u_{\boldsymbol{q}},v_{\boldsymbol{q}},w_{\boldsymbol{q}}$ and
$x_{\boldsymbol{q}}$ can be defined as real functions, thus the previous
relations for the mean-field coefficients are simplified in the following
form,

$\displaystyle{\cal C}_{\alpha\alpha,{\boldsymbol{q}}}$
$\displaystyle=\Delta^{-}_{\boldsymbol{q}}|u_{\boldsymbol{q}}|^{2}+\Delta^{-}_{-\boldsymbol{q}}|x_{\boldsymbol{q}}|^{2}+2\,\text{Re}\left[u_{\boldsymbol{q}}x_{\boldsymbol{q}}\,\big{(}\xi-\text{Re}\,\Gamma_{q}\big{)}\right],$
(30a) $\displaystyle{\cal C}_{\beta\beta,{\boldsymbol{q}}}$
$\displaystyle=\Delta^{+}_{\boldsymbol{q}}|v_{\boldsymbol{q}}|^{2}+\Delta^{+}_{\boldsymbol{q}}|w_{\boldsymbol{q}}|^{2}+2\,\text{Re}\left[v_{\boldsymbol{q}}w_{\boldsymbol{q}}\,\big{(}\xi+\text{Re}\,\Gamma_{q}\big{)}\right],$
$\displaystyle{\cal C}^{\prime}_{\alpha\alpha,{\boldsymbol{q}}}$
$\displaystyle=2\Delta^{-}_{\boldsymbol{q}}u_{\boldsymbol{q}}{x}_{\boldsymbol{q}}+\Lambda^{-}_{\boldsymbol{q}}u^{2}_{\boldsymbol{q}}+\bar{\Lambda}^{-}_{\boldsymbol{q}}{x}^{2}_{\boldsymbol{q}},$
$\displaystyle{\cal C}^{\prime}_{\beta\beta,{\boldsymbol{q}}}$
$\displaystyle=2\Delta^{+}_{\boldsymbol{q}}v_{\boldsymbol{q}}{w}_{\boldsymbol{q}}+\Lambda^{+}_{\boldsymbol{q}}v^{2}_{\boldsymbol{q}}+\bar{\Lambda}^{+}_{\boldsymbol{q}}{w}^{2}_{\boldsymbol{q}},$
$\displaystyle{\cal D}_{\alpha\beta,{\boldsymbol{q}}}$
$\displaystyle=\Gamma_{\boldsymbol{q}}x_{\boldsymbol{q}}v_{\boldsymbol{q}}-\bar{\Gamma}_{\boldsymbol{q}}{u}_{\boldsymbol{q}}{w}_{\boldsymbol{q}}+\Theta_{-\boldsymbol{q}}x_{\boldsymbol{q}}{w}_{\boldsymbol{q}}-\Theta_{\boldsymbol{q}}v_{\boldsymbol{q}}{u}_{\boldsymbol{q}}+2\,\xi
v_{\boldsymbol{q}}\,x_{\boldsymbol{q}}+2\,\bar{\xi}\,u_{\boldsymbol{q}}\,w_{\boldsymbol{q}},$
$\displaystyle{\cal D}^{\prime}_{\alpha\beta,{\boldsymbol{q}}}$
$\displaystyle=\bar{\Gamma}_{\boldsymbol{q}}{x}_{\boldsymbol{q}}{w}_{\boldsymbol{q}}-\Gamma_{\boldsymbol{q}}u_{\boldsymbol{q}}v_{\boldsymbol{q}}-{\Theta}_{-\boldsymbol{q}}{x}_{\boldsymbol{q}}v_{\boldsymbol{q}}+\Theta_{\boldsymbol{q}}{w}_{\boldsymbol{q}}u_{\boldsymbol{q}}+2\xi^{\prime}u_{\boldsymbol{q}}v_{\boldsymbol{q}}+2\,\bar{\xi}^{\prime}\,x_{\boldsymbol{q}}w_{\boldsymbol{q}}.$

where we define,

$\displaystyle\Delta^{\pm}_{\boldsymbol{q}}$
$\displaystyle=A_{a}\pm\text{Re}\left[A^{\prime}_{ab,{\boldsymbol{q}}}\right],$
(31a) $\displaystyle\xi$ $\displaystyle=A_{a}^{\prime}+A_{b}^{\prime},\hskip
85.35826pt\xi^{\prime}=A_{a}^{\prime}-A_{b}^{\prime},$ (31b)
$\displaystyle\Gamma_{\boldsymbol{q}}$
$\displaystyle=(A_{ab,-{\boldsymbol{q}}}-A_{ab,{\boldsymbol{q}}}),\hskip
56.9055pt\Gamma_{\boldsymbol{q}}^{\prime}=(A_{ab,-{\boldsymbol{q}}}+A_{ab,{\boldsymbol{q}}}),$
(31c) $\displaystyle\Lambda^{\pm}_{\boldsymbol{q}}$ $\displaystyle=\xi\pm
A_{ab,{\boldsymbol{q}}},$ (31d) $\displaystyle\Theta_{\boldsymbol{q}}$
$\displaystyle=(A^{\prime}_{ab,{\boldsymbol{q}}}-\bar{A}^{\prime}_{ab,{\boldsymbol{q}}}).$
(31e)

###  B.3 Eigenenergies of a biaxial AF system

The eigenenergy of a biaxial AF Hamiltonian in the presence of nonlinear
interaction corrections, modelled by an effective mean-field Hamiltonian
$H=H^{(2)}_{\text{m}}+H^{(4)}_{\text{MF}}$, where $H=H^{(2)}_{\text{m}}$ is
given by Eq. (3) and $H^{(4)}_{\text{MF}}$ is given by Eq. (10); are given by,

$\displaystyle 2\tilde{\epsilon}^{2}_{\mu,\nu}={\cal X}$
$\displaystyle\pm\left[{{\cal X}^{2}-\left({\cal D}^{\prime
2}_{\alpha\beta}\left(0.25\bar{{\cal D}}^{\prime 2}_{\alpha\beta}-\bar{\cal
C}^{\prime}_{\alpha\alpha}\bar{\cal
C}^{\prime}_{\beta\beta}\right)+4\left(\epsilon_{\alpha}+0.5{\cal
C}_{\alpha\alpha}\right)^{2}\left(\left(\epsilon_{\beta}+0.5{\cal
C}_{\beta\beta}\right)^{2}-{\cal C}^{\prime}_{\beta\beta}\bar{\cal
C}^{\prime}_{\beta\beta}\right)\right)}\right.$ $\displaystyle\left.+{\cal
C}^{\prime}_{\alpha\alpha}\left(4\left(\epsilon_{\beta}+0.5{\cal
C}_{\beta\beta}\right)^{2}\bar{\cal C}^{\prime}_{\alpha\alpha}+\bar{{\cal
D}}_{\alpha\beta}^{\prime 2}{\cal C}^{\prime}_{\beta\beta}-4\bar{\cal
C}^{\prime}_{\alpha\alpha}|{\cal
C}^{\prime}_{\beta\beta}|^{2}-2\left(\epsilon_{\beta}+0.5{\cal
C}_{\beta\beta}\right)\bar{{\cal D}}^{\prime}_{\alpha\beta}{\cal
D}_{\alpha\beta}+\bar{\cal C}^{\prime}_{\beta\beta}{\cal
D}_{\alpha\beta}^{2}\right)\right.$ $\displaystyle\left.-2{\cal
D}^{\prime}_{\alpha\beta}\bar{\cal
D}_{\alpha\beta}\left(\left(\epsilon_{\beta}+0.5{\cal
C}_{\beta\beta}\right)\bar{\cal C}^{\prime}_{\alpha\alpha}-0.25\bar{{\cal
D}}^{\prime}_{\alpha\beta}{\cal D}_{\alpha\beta}\right)-\left(-\bar{\cal
C}^{\prime}_{\alpha\alpha}{\cal C}^{\prime}_{\beta\beta}+0.25{\cal
D}^{2}_{\alpha\beta}\right)\bar{\cal D}^{2}_{\alpha\beta}\right.$
$\displaystyle\left.-2\left(\epsilon_{\alpha}+0.5{\cal
C}_{\alpha\alpha}\right)\left({\cal D}^{\prime}_{\alpha\beta}\bar{\cal
C}^{\prime}_{\beta\beta}{\cal D}_{\alpha\beta}+\bar{{\cal
D}}^{\prime}_{\alpha\beta}{\cal C}^{\prime}_{\beta\beta}\bar{\cal
D}_{\alpha\beta}-\left(\epsilon_{\beta}+0.5{\cal
C}_{\beta\beta}\right)\left(|{\cal D}^{\prime}_{\alpha\beta}|^{2}+|{\cal
D}_{\alpha\beta}|^{2}\right)\right)\right]^{\frac{1}{2}},$ (32)

with ${\cal X}=\left(\epsilon_{\alpha}+0.5{\cal
C}_{\alpha\alpha}\right)^{2}+\left(\epsilon_{\beta}+0.5{\cal
C}_{\beta\beta}\right)^{2}+0.5\left(\left|{\cal
D}_{\alpha\beta}\right|^{2}-\left|{\cal
D}^{\prime}_{\alpha\beta}\right|^{2}\right)-\left|{\cal
C}^{\prime}_{\alpha\alpha}\right|^{2}-\left|{\cal
C}^{\prime}_{\beta\beta}\right|^{2}$. In square lattices, as we have already
discussed in the main text, the interband scattering processes are absent
${\cal D}_{\alpha\beta}={\cal D}^{\prime}_{\alpha\beta}=0$, and thus Eq. (B.3)
reduces to Eq. (17) in the main text. The result of Eq. (B.3), is computed for
cubic lattices and presented at Fig. 6 for different temperatures.

Figure 6: Dispersion relation along the high-symmetry lines for temperature-
dependent magnonic eigenmodes of a 3D cubic biaxial AF insulator. The
eigenenergies for each mode, $\epsilon_{\mu}$ (left panel) and
$\epsilon_{\nu}$ (right panel), are presented for various temperatures, using
Eq. (B.3).