Thermal squeezing and nonlinear spectral shift of magnons in antiferromagnetic insulators

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 nanotechnologies [1–5]. 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 polarization degree of freedom of magnons [6–11]. Of particular interest are the electrical insulator class of AF materials [10, 12–19], 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 shown [3, 19, 20].

Magnonics is an emerging field in spintronics and magnetic-based nanotechnology [12] that exhibits various features from quantum coherent phenomena [21–24], chiral magnonic spin transport, magnonic crystals [25–28], and ultrafast magnetization dynamics [29]. 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 [30, 31] through magnon relaxation and spin conductivity, as well as close to magnetic phase transitions. Moreover, nonlinear magnon interactions are essential to stabilize magnetic droplets [32] and Bose–Einstein condensation of magnons [23, 33–35]. It is also an important ingredient in emergent magnonic topological phases [36]. A recent study demonstrated that nonlinear magnon interactions of Dirac magnons in honeycomb ferromagnets strongly renormalize the magnon band structure and magnon lifetimes [37]. In AF systems, nonlinear magnon interactions also have significant effects even at low temperatures due to zero-point fluctuations [38]. Additionally, two-magnon spectra of AF systems are strongly influenced by magnon–magnon interactions [39].

The recent discovery of two-dimensional (2D) materials establishes a new platform for the study of nonlinear magnon interactions in different lattice structures [31, 40–43]. 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 [44, 45]. Therefore, the effect of nonlinear 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 nonlinear magnon interactions in AF systems within a mean-field approach. In references [46, 47], 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 reference [48] 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. For a recent review on AF systems look at reference [20].

Nonlinear magnon interactions, even in simple AF systems, are very complicated. To the best of our knowledge, a systematic investigation of nonlinear interactions on spectrum and polarization of thermal magnons in AF materials remains absent [49]. In AF systems, the underlying crystal lattice structure affects the magnetic ground states and interactions, and therefore, one needs to study each case separately [50].

In this paper, we first develop a formalism to compute the nonlinear interactions of thermal magnons in both uniaxial easy-axis and biaxial easy-plane AF insulators in arbitrary dimensions using a self-consistent Hartree–Fock mean-field theory. Next, we investigate the effect of thermal nonlinear interactions on bare magnon dispersion and AF resonance frequencies, which are directly accessible experimentally [51]. Later, we investigate the effect of nonlinear interactions on the polarization of magnon modes. For numerical calculations, we consider collinear AF materials with easy-axis and easy-plane 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 section 2, we introduce the spin Hamiltonian and review the linear spin-wave theory. In section 3, we present a Hartree–Fock mean-field theory for thermal magnon interactions and compute mean-field coefficients. In section 4, we compute the nonlinear spectral shift of magnons in easy-axis and easy-plane AF systems and present numerical calculations for two specific lattice models. In section 5, we investigate the effect of nonlinear thermal interactions on the polarization of magnons. We conclude in section 6 with a discussion of our results. In the appendices, we present the technical details of the calculations.

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, [6, 20],

$$\begin{equation}{H}_{s}=J\sum\limits _{\langle ij\rangle }{\mathbf{S}}_{i}\cdot {\mathbf{S}}_{j}+\sum\limits _{i}{\mathcal{K}}_{x}{({S}_{i}^{x})}^{2}-\sum\limits _{i}{\mathcal{K}}_{z}{({S}_{i}^{z})}^{2},\end{equation} \tag{ 1 }$$

where ⟨ij⟩ represents summation over nearest-neighbor sites, i and j, and we consider the local spins have the same length S at each sublattice $\vert {\mathbf{S}}_{i\in \left\{\mathcal{A}\right\}}\vert =\vert {\mathbf{S}}_{i\in \left\{\mathcal{B}\right\}}\vert =S$. The Hamiltonian includes the isotropic quantum AF Heisenberg exchange interaction, J > 0, and single-ion 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 easy axis case, i.e., ${\mathcal{K}}_{x}=0$, both the spin Hamiltonian and the ground state have U(1) symmetry while in the biaxial easy-plane 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, equation (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 and the eigenmodes are not degenerate anymore [52].

To study spin excitations, we write the spin Hamiltonian, equation (1), in terms of so-called deviation bosonic operators using the Holstein–Primakoff (HP) transformation [53]. We assume the ground state of the AF spin Hamiltonian is a collinear Néel state with ${\mathbf{S}}_{i\in \left\{\mathcal{A}\right\}}=(0,0,S)$ and ${\mathbf{S}}_{j\in \left\{\mathcal{B}\right\}}=(0,0,-S)$; and thus the HP transformation in the two-sublattice AF system reads,

$$\begin{equation}{S}_{i}^{z}=S-{a}_{i}^{{\dagger}}{a}_{i},\quad {S}_{i}^{-}={a}_{i}^{{\dagger}}\sqrt{2S-{a}_{i}^{{\dagger}}{a}_{i}},\end{equation} \tag{ 2a }$$

$$\begin{equation}{S}_{j}^{z}=-S+{b}_{j}^{{\dagger}}{b}_{j},\quad {S}_{j}^{-}=\sqrt{2S-{b}_{j}^{{\dagger}}{b}_{j}}\,{b}_{j},\end{equation} \tag{ 2b }$$

and ${S}_{i(j)}^{+}={\left({S}_{i(j)}^{-}\right)}^{{\dagger}}$, where i and j denote sites in the sublattices $\mathcal{A}$ and $\mathcal{B}$, respectively. The deviation operators a (b) and a^† (b^†) 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 [54]. Inserting equations (2a) and (2b) into the AF spin Hamiltonian, equation (1), and under the assumption of small deviation number $\langle {a}_{i}^{{\dagger}}{a}_{i}\rangle ,\langle {b}_{i}^{{\dagger}}{b}_{i}\rangle \ll 2S$, with ⟨⋅⟩ 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,

$$\begin{equation}{H}_{b}={E}_{0}^{c}+{H}_{b}^{(2)}+{H}_{b}^{(4)}+\cdots \,,\end{equation} \tag{ 3 }$$

where ellipsis stands for higher order interactions. The classical AF ground-state energy is given by,

$$\begin{equation}{E}_{0}^{c}=NS\left[\left.-(zJ/2+{\mathcal{K}}_{z})S+{\mathcal{K}}_{x}/2\right)\right],\end{equation} \tag{ 4 }$$

where z denotes the coordination number and N is the total number of sites. ${H}_{b}^{(2)}$ and ${H}_{b}^{(4)}$ represent the noninteracting (quadric) and interacting (quartic) bosonic Hamiltonians, respectively.

The noninteracting Hamiltonian, consists of quadratic bosonic terms, in the momentum space, reads,

$$\begin{align}{H}_{b}^{(2)}& =S\sum\limits _{\boldsymbol{q}}\left[\right.(zJ+{\mathcal{K}}_{x}+2{\mathcal{K}}_{z})({a}_{\boldsymbol{q}}^{{\dagger}}{a}_{\boldsymbol{q}}+{b}_{\boldsymbol{q}}^{{\dagger}}{b}_{\boldsymbol{q}})+zJ{\gamma }_{-\boldsymbol{q}}{a}_{\boldsymbol{q}}{b}_{-\boldsymbol{q}}\\ & \quad +zJ{\gamma }_{\boldsymbol{q}}{b}_{-\boldsymbol{q}}^{{\dagger}}{a}_{\boldsymbol{q}}^{{\dagger}}+\frac{{\mathcal{K}}_{x}}{2}\left({a}_{\boldsymbol{q}}{a}_{-\boldsymbol{q}}+{b}_{\boldsymbol{q}}{b}_{-\boldsymbol{q}}+{a}_{\boldsymbol{q}}^{{\dagger}}{a}_{-\boldsymbol{q}}^{{\dagger}}+{b}_{\boldsymbol{q}}^{{\dagger}}{b}_{-\boldsymbol{q}}^{{\dagger}}\right)\left.\right],\end{align} \tag{ 5 }$$

where ${a}_{i}=\sqrt{2/N}{\sum }_{\boldsymbol{q}}{\text{e}}^{\text{i}\boldsymbol{q}\cdot {\mathbf{r}}_{i}}{a}_{\boldsymbol{q}}$ and ${b}_{j}=\sqrt{2/N}{\sum }_{\boldsymbol{q}}{\text{e}}^{\text{i}\boldsymbol{q}\cdot {\mathbf{r}}_{j}}{b}_{\boldsymbol{q}}$, with $\left\{i,j\right\}\in \left\{\mathcal{A},\mathcal{B}\right\}$. The lowest order interacting Hamiltonian, in our model, consists of quartic bosonic terms and describes the four-magnon scattering processes ${H}_{b}^{(4)}$. In section 3, we will investigate these nonlinear terms in detail. In general, in AF spin systems, three-magnon interactions, consists of cubic bosonic terms, can be appeared in the presence of some spin interactions, like DMI and dipole–dipole interactions, as well as magnetic frustrations. These interactions break the conservation of total magnon numbers. However, in our model, equation (1), due to the absence of such spin interactions, three-magnon interactions are forbidden [38].

2.1. Linear spin-wave theory

To investigate the effect of nonlinear four-magnon interactions on magnon dispersion and polarization, we should first find the band structure and eigenstates of the noninteracting magnon modes. At a magnon density, we can disregard nonlinear magnon interactions. The quadratic bosonic Hamiltonian ${H}_{b}^{(2)}$ can be diagonalized via a standard canonical Bogolioubov transformation that maps the boson operators on two AF sublattices {a _q , b _q }, into the new bosonic magnon eigenmodes {α _q , β _q } where q is the magnon wavevector. In the rest of this section, we find the eigensystems of uniaxial and biaxial AF systems, separately.

2.1.1. Uniaxial easy-axis AF system

To diagonalize the noninteracting bosonic Hamiltonian (5) in the uniaxial easy-axis limit $({\mathcal{K}}_{x}=0)$, we perform a canonical transformation using the following 2 × 2 Bogolioubov transformation, [20, 55],

$$\begin{equation}\left(\begin{matrix}\hfill {a}_{\boldsymbol{q}}\hfill \\ \hfill {b}_{-\boldsymbol{q}}^{{\dagger}}\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill {\bar{u}}_{\boldsymbol{q}}^{\mathrm{u}}\hfill & \hfill -{\bar{v}}_{\boldsymbol{q}}^{\mathrm{u}}\hfill \\ \hfill -{v}_{\boldsymbol{q}}^{\mathrm{u}}\hfill & \hfill {u}_{\boldsymbol{q}}^{\mathrm{u}}\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {\alpha }_{\boldsymbol{q}}^{\mathrm{u}}\hfill \\ \hfill {\beta }_{\boldsymbol{q}}^{\mathrm{u}{\dagger}}\hfill \end{matrix}\right),\end{equation} \tag{ 6 }$$

where $\vert {u}_{\boldsymbol{q}}^{\mathrm{u}}{\vert }^{2}-\vert {v}_{\boldsymbol{q}}^{\mathrm{u}}{\vert }^{2}=1$. We define ${u}_{\boldsymbol{q}}^{\mathrm{u}}=\mathrm{cosh}({\theta }_{\boldsymbol{q}}/2)$ and ${v}_{\boldsymbol{q}}^{\mathrm{u}}=\mathrm{sinh}({\theta }_{\boldsymbol{q}}/2)$. Inserting equation (6) into the Hamiltonian (5), we obtain the following Bogoliubov coefficients for the uniaxial easy-axis system,

$$\begin{equation}{u}_{\boldsymbol{q}}^{\mathrm{u}}=\sqrt{\frac{zJS+2{\mathcal{K}}_{z}S+{{\epsilon}}_{\alpha ,\boldsymbol{q}}^{u}}{2{{\epsilon}}_{\alpha ,\boldsymbol{q}}^{u}}},\end{equation} \tag{ 7a }$$

$$\begin{equation}{v}_{\boldsymbol{q}}^{\mathrm{u}}=\sqrt{\frac{zJS+2{\mathcal{K}}_{z}S-{{\epsilon}}_{\beta ,\boldsymbol{q}}^{u}}{2{{\epsilon}}_{\beta ,\boldsymbol{q}}^{u}}}.\end{equation} \tag{ 7b }$$

The degenerate eigenenergies of the circularly polarized eigenmodes χ ∈ {α, β} of easy-axis AF system are,

$$\begin{equation}{{\epsilon}}_{\chi ,\boldsymbol{q}}^{\mathrm{u}}=S\sqrt{{\left(zJ+2{\mathcal{K}}_{z}\right)}^{2}-{\left(zJ\vert {\gamma }_{\boldsymbol{q}}\vert \right)}^{2}},\end{equation} \tag{ 8 }$$

where ${\gamma }_{\boldsymbol{q}}={z}^{-1}{\sum }_{i=1}^{z}{\text{e}}^{\text{i}\boldsymbol{q}\cdot {\boldsymbol{\delta }}_{i}}$ is the structure factor and δ _i denotes the nearest-neighbor vectors.

2.1.2. Biaxial easy-plane AF system

In the case of biaxial AF systems, we use the following 4 × 4 Bogolioubov transformations [20, 56],

$$\begin{equation}\left(\begin{matrix}\hfill {a}_{\boldsymbol{q}}\hfill \\ \hfill {b}_{\boldsymbol{q}}\hfill \\ \hfill {a}_{-\boldsymbol{q}}^{{\dagger}}\hfill \\ \hfill {b}_{-\boldsymbol{q}}^{{\dagger}}\hfill \end{matrix}\right)=\left(\begin{matrix}\hfill {u}_{\boldsymbol{q}}\hfill & \hfill \hfill & \hfill {v}_{\boldsymbol{q}}\hfill & \hfill \hfill & \hfill {x}_{\boldsymbol{q}}\hfill & \hfill \hfill & \hfill {w}_{\boldsymbol{q}}\hfill \\ \hfill -{u}_{\boldsymbol{q}}\hfill & \hfill \hfill & \hfill {v}_{\boldsymbol{q}}\hfill & \hfill \hfill & \hfill -{x}_{\boldsymbol{q}}\hfill & \hfill \hfill & \hfill {w}_{\boldsymbol{q}}\hfill \\ \hfill {\bar{x}}_{\boldsymbol{q}}\hfill & \hfill \hfill & \hfill {\bar{w}}_{\boldsymbol{q}}\hfill & \hfill \hfill & \hfill {\bar{u}}_{\boldsymbol{q}}\hfill & \hfill \hfill & \hfill {\bar{v}}_{\boldsymbol{q}}\hfill \\ \hfill -{\bar{x}}_{\boldsymbol{q}}\hfill & \hfill \hfill & \hfill {\bar{w}}_{\boldsymbol{q}}\hfill & \hfill \hfill & \hfill -{\bar{u}}_{\boldsymbol{q}}\hfill & \hfill \hfill & \hfill {\bar{v}}_{\boldsymbol{q}}\hfill \end{matrix}\right)\left(\begin{matrix}\hfill {\alpha }_{\boldsymbol{q}}\hfill \\ \hfill {\beta }_{\boldsymbol{q}}\hfill \\ \hfill {\alpha }_{-\boldsymbol{q}}^{{\dagger}}\hfill \\ \hfill {\beta }_{-\boldsymbol{q}}^{{\dagger}}\hfill \end{matrix}\right),\end{equation} \tag{ 9 }$$

where the Bogoliubov coefficients are given by,

$$\begin{equation}{u}_{\boldsymbol{q}}=\frac{1}{2}\sqrt{\frac{zJS+{\mathcal{K}}_{x}S+2{\mathcal{K}}_{z}S+{{\epsilon}}_{\boldsymbol{q}}^{\alpha }}{{{\epsilon}}_{\boldsymbol{q}}^{\alpha }}},\end{equation} \tag{ 10a }$$

$$\begin{equation}{v}_{\boldsymbol{q}}=-\frac{1}{2}\sqrt{\frac{zJS+{\mathcal{K}}_{x}S+2{\mathcal{K}}_{z}S+{{\epsilon}}_{\boldsymbol{q}}^{\beta }}{{{\epsilon}}_{\boldsymbol{q}}^{\beta }}},\end{equation} \tag{ 10b }$$

$$\begin{equation}{x}_{\boldsymbol{q}}=\frac{1}{2}\sqrt{\frac{zJS+{\mathcal{K}}_{x}S+2{\mathcal{K}}_{z}S-{{\epsilon}}_{\boldsymbol{q}}^{\alpha }}{{{\epsilon}}_{\boldsymbol{q}}^{\alpha }}},\end{equation} \tag{ 10c }$$

$$\begin{equation}{w}_{\boldsymbol{q}}=\frac{1}{2}\sqrt{\frac{zJS+{\mathcal{K}}_{x}S+2{\mathcal{K}}_{z}S-{{\epsilon}}_{\boldsymbol{q}}^{\beta }}{{{\epsilon}}_{\boldsymbol{q}}^{\beta }}}.\end{equation} \tag{ 10d }$$

The eigenenergies of two magnon modes χ ∈ {α, β} in easy-plane AF system are,

$$\begin{equation}{{\epsilon}}_{\boldsymbol{q}}^{\chi }=S\sqrt{{\left(zJ+{\mathcal{K}}_{x}+2{\mathcal{K}}_{z}\right)}^{2}-{\left(zJ\vert {\gamma }_{\boldsymbol{q}}\vert \mp {\mathcal{K}}_{x}\right)}^{2}},\end{equation} \tag{ 11 }$$

where sign −(+) refers to α(β) magnon mode. The eignenenergies of the biaxial easy-plane case reduces to the uniaxial easy-axis limit when ${\mathcal{K}}_{x}=0$, as it should be.

In the 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 [7, 20]. The presence of a hard-axis anisotropy $({\mathcal{K}}_{x}\ne 0)$ removes the degeneracy of magnon eigenstates and causes the two eigenmodes become elliptically polarized [7, 20]. The latter is a consequence of the breaking of U(1) symmetry in the ground state of the spin Hamiltonian 1. In the presence of the long-range dipolar interactions that also breaks U(1) symmetry, the magnon eigenmodes of a uniaxial easy-axis are elliptically polarized [57]. In the present study, we ignore the effect of dipolar interaction since they are negligible in most of AF systems.

After diagonalization, the total noninteracting AF Hamiltonian in the new eigenbasis is given by,

$$\begin{equation}{H}_{m}={E}_{0}+{H}_{m}^{(2)},\end{equation} \tag{ 12 }$$

$$\begin{equation}{E}_{0}={E}_{0}^{c}+{E}_{0}^{q},\end{equation} \tag{ 13 }$$

where E₀ is the total AF ground-state energy, and the quantum correction to the classic AF ground-state energy, ${E}_{0}^{c}$, reads,

$$\begin{equation}{E}_{0}^{q}=S\sum\limits _{\mathbf{q}}\left\{2(zJ+{\mathcal{K}}_{x}+2{\mathcal{K}}_{z})(\vert {x}_{\mathbf{q}}{\vert }^{2}+\vert {w}_{\mathbf{q}}{\vert }^{2})+2\,\text{Re}\left[zJ{\gamma }_{-\mathbf{q}}(-{u}_{\mathbf{q}}{x}_{\mathbf{q}}+{v}_{\mathbf{q}}{w}_{\mathbf{q}})+{\mathcal{K}}_{x}({u}_{\mathbf{q}}{x}_{\mathbf{q}}+{v}_{\mathbf{q}}{w}_{\mathbf{q}})\right]\right\}.\end{equation} \tag{ 14 }$$

This quantum zero-point energy 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 [58]. The quadratic part of the AF Hamiltonian that describes the dynamics of magnons is,

$$\begin{equation}{H}_{m}^{(2)}=\sum\limits _{\boldsymbol{q}}\left[{{\epsilon}}_{\boldsymbol{q}}^{\alpha }{\alpha }_{\boldsymbol{q}}^{{\dagger}}{\alpha }_{\boldsymbol{q}}+{{\epsilon}}_{\boldsymbol{q}}^{\beta }{\beta }_{\boldsymbol{q}}^{{\dagger}}{\beta }_{\boldsymbol{q}}\right],\end{equation} \tag{ 15 }$$

where the bare eigenenergy for the magnon mode χ ∈ {α, β} is given by equation (11).

The structure factors for 2D square (z = 4) and hexagonal (z = 3) lattices are, ${\gamma }_{\boldsymbol{q}}=2\left(\right.\cos ({q}_{x}/{q}_{m})+\cos ({q}_{y}/{q}_{m})\left.\right)/z$ and ${\gamma }_{\boldsymbol{q}}={\text{e}}^{\text{i}{q}_{x}/{q}_{m}}\left(1+2\,{\text{e}}^{-\text{i}3{q}_{x}/2{q}_{m}}\,\cos (\sqrt{3}{q}_{y}/2{q}_{m})\right)/z$, respectively, where q_m = 1/a and a is the lattice constant. In the long-wavelength limit q → 0, the structure factors of both hexagonal and square lattices reduce to ${\gamma }_{\boldsymbol{q}}\simeq 1-{q}^{2}/4{q}_{m}^{2}$; and thus the mode-dependent eigenenergies becomes to ${{\epsilon}}_{\mathbf{q}}^{\alpha }=S\sqrt{4({\mathcal{K}}_{z}+{\mathcal{K}}_{x})(zJ+{\mathcal{K}}_{z})+zJ{q}^{2}(zJ-{\mathcal{K}}_{x})/2{q}_{m}^{2}}$ and ${{\epsilon}}_{\mathbf{q}}^{\beta }=S\sqrt{4{\mathcal{K}}_{z}(zJ+{\mathcal{K}}_{z}+{\mathcal{K}}_{x})+zJ{q}^{2}(zJ+{\mathcal{K}}_{x})/2{q}_{m}^{2}}$. In other words, around the center of the Brillouin zone, called Γ-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}}_{\mathbf{q}=0}^{\alpha }$ and ${{\epsilon}}_{\mathbf{q}=0}^{\beta }$ 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.

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,

$$\begin{equation}{H}_{b}^{(4)}=\frac{-J}{4}\sum\limits _{\langle i,j\rangle }\left[{a}_{i}^{{\dagger}}{a}_{i}{b}_{j}^{{\dagger}}{b}_{j}+{a}_{i}{b}_{j}^{{\dagger}}{b}_{j}{b}_{j}+(a{\leftrightarrow}b)\right]-\frac{{\mathcal{K}}_{x}}{8}\sum\limits _{i}\left[{\Delta}\mathcal{K}{a}_{i}^{{\dagger}}{a}_{i}^{{\dagger}}{a}_{i}{a}_{i}+{a}_{i}^{{\dagger}}{a}_{i}{a}_{i}{a}_{i}+(a{\leftrightarrow}b)\right]+\text{H.c.},\end{equation} \tag{ 16 }$$

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 equation (16), 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 [38, 59].

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 [60, 61] to reduce equation (16) into an effective noninteracting Hamiltonian in terms of quadratic magnon operators. In the mean-field treatment, we start out by substituting ${\xi }_{i}^{{\dagger}}{\xi }_{j}^{\prime }={\langle {\xi }_{i}^{{\dagger}}{\xi }_{j}^{\prime }\rangle }_{\text{th}}+{{\Xi}}_{\xi {\xi }^{\prime }}$ into the Hamiltonian (16), where ${\xi }_{i},{\xi }_{i}^{\prime }\in \left\{{a}_{i},{b}_{i},{a}_{i}^{{\dagger}},{b}_{i}^{{\dagger}}\right\}$ and ⟨⋅⟩_th is a thermal average. We expand up to the first order in the field ${{\Xi}}_{\xi {\xi }^{\prime }}\equiv {\xi }_{i}^{{\dagger}}{\xi }_{j}^{\prime }-{\langle {\xi }_{i}^{{\dagger}}{\xi }_{j}^{\prime }\rangle }_{\text{th}}$, and express the operators in the basis of α- and β-magnons through the Bogoliubov transformation. In the leading-order correction, it is assumed that mean-field parameter Ξ_ξξ' 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 [62]. In the Hartree–Fock approximation, we only keep the diagonal terms and assume for each bosonic operator χ _q ∈ {α _q , β _q }, so that ${\langle {\chi }_{\boldsymbol{q}}^{{\dagger}}{\chi }_{{\boldsymbol{q}}^{\prime }}\rangle }_{\text{th}}={\delta }_{\boldsymbol{q},{\boldsymbol{q}}^{\prime }}{n}_{\chi ,\boldsymbol{q}}$[42, 43]. The thermodynamic equilibrium of the AF magnon number at each magnon mode χ ∈ {α, β} is defined through the Bose–Einstein distribution function ${n}_{\chi ,\boldsymbol{q}}={({\mathrm{e}}^{{{\epsilon}}_{\chi ,\boldsymbol{q}}/{k}_{\mathrm{B}}T}-1)}^{-1}$, with the Boltzmann constant k_B, and the temperature T. The chemical potential of magnons is zero at thermal equilibrium [63].

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

$$\begin{equation}{\mathbb{H}}^{(4)}=\left(\begin{matrix}\hfill {\mathcal{C}}_{\alpha \alpha }\hfill & \hfill {\mathcal{D}}_{\alpha \beta }\hfill & \hfill 2{\bar{\mathcal{C}}}_{\alpha \alpha }^{\prime }\hfill & \hfill {\bar{\mathcal{D}}}_{\alpha \beta }^{\prime }\hfill \\ \hfill {\bar{\mathcal{D}}}_{\alpha \beta }\hfill & \hfill {\mathcal{C}}_{\beta \beta }\hfill & \hfill {\bar{\mathcal{D}}}_{\alpha \beta }^{\prime }\hfill & \hfill 2{\bar{\mathcal{C}}}_{\beta \beta }^{\prime }\hfill \\ \hfill 2{\mathcal{C}}_{\alpha \alpha }^{\prime }\hfill & \hfill {\mathcal{D}}_{\alpha \beta }^{\prime }\hfill & \hfill {\mathcal{C}}_{\alpha \alpha }\hfill & \hfill {\bar{\mathcal{D}}}_{\alpha \beta }\hfill \\ \hfill {\mathcal{D}}_{\alpha \beta }^{\prime }\hfill & \hfill 2{\mathcal{C}}_{\beta \beta }^{\prime }\hfill & \hfill {\mathcal{D}}_{\alpha \beta }\hfill & \hfill {\mathcal{C}}_{\beta \beta }\hfill \end{matrix}\right),\end{equation} \tag{ 17 }$$

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 A2. 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 ${\mathcal{C}}_{\chi \chi }$ and ${\mathcal{C}}_{\chi \chi }^{\prime }$, describing intraband processes of magnons at the χ-band, are given by equation (A7) in appendix A2.

$$\begin{equation}{\mathcal{C}}_{\alpha \alpha ,\boldsymbol{q}}={{\Delta}}_{\boldsymbol{q}}^{-}\vert {u}_{\boldsymbol{q}}{\vert }^{2}+{{\Delta}}_{-\boldsymbol{q}}^{-}\vert {x}_{\boldsymbol{q}}{\vert }^{2}+2\,\text{Re}\left[{u}_{\boldsymbol{q}}{x}_{\boldsymbol{q}}\,\left(\xi -\text{Re}\,{{\Gamma}}_{q}\right)\right],\end{equation} \tag{ 18a }$$

$$\begin{equation}{\mathcal{C}}_{\beta \beta ,\boldsymbol{q}}={{\Delta}}_{\boldsymbol{q}}^{+}\vert {v}_{\boldsymbol{q}}{\vert }^{2}+{{\Delta}}_{\boldsymbol{q}}^{+}\vert {w}_{\boldsymbol{q}}{\vert }^{2}+2\,\text{Re}\left[{v}_{\boldsymbol{q}}{w}_{\boldsymbol{q}}\,\left(\xi +\text{Re}\,{{\Gamma}}_{q}\right)\right],\end{equation} \tag{ 18b }$$

$$\begin{equation}{\mathcal{C}}_{\alpha \alpha ,\boldsymbol{q}}^{\prime }=2{{\Delta}}_{\boldsymbol{q}}^{-}{u}_{\boldsymbol{q}}{x}_{\boldsymbol{q}}+{{\Lambda}}_{\boldsymbol{q}}^{-}{u}_{\boldsymbol{q}}^{2}+{\bar{{\Lambda}}}_{\boldsymbol{q}}^{-}{x}_{\boldsymbol{q}}^{2},\end{equation} \tag{ 18c }$$

$$\begin{equation}{\mathcal{C}}_{\beta \beta ,\boldsymbol{q}}^{\prime }=2{{\Delta}}_{\boldsymbol{q}}^{+}{v}_{\boldsymbol{q}}{w}_{\boldsymbol{q}}+{{\Lambda}}_{\boldsymbol{q}}^{+}{v}_{\boldsymbol{q}}^{2}+{\bar{{\Lambda}}}_{\boldsymbol{q}}^{+}{w}_{\boldsymbol{q}}^{2},\end{equation} \tag{ 18d }$$

where ${{\Delta}}_{\boldsymbol{q}}^{\pm }$, ξ, Γ _q and ${{\Lambda}}_{\boldsymbol{q}}^{\pm }$ are defined in equation (A8). These mean-field coefficients parameterize the amplitude of the intra- and inter-band magnon interactions. The relevant Bogoliubov coefficients u _q , v _q , x _q and w _q are defined by equations (10a)–(10d). Interband scattering of magnons are characterized by the components ${\mathcal{D}}_{\alpha \beta }$ and ${\mathcal{D}}_{\alpha \beta }^{\prime }$, and read,

$$\begin{equation}{\mathcal{D}}_{\alpha \beta ,\boldsymbol{q}}={{\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}},\end{equation} \tag{ 19a }$$

$$\begin{equation}{\mathcal{D}}_{\alpha \beta ,\boldsymbol{q}}^{\prime }={\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}},\end{equation} \tag{ 19b }$$

where ${{\Gamma}}_{\boldsymbol{q}}^{\prime }$, Θ _q and ξ' can be found in equation (A8). Through the interband process, α- and β-magnon modes become coupled and thus, hybridized. The definition of ${A}_{a},\,{A}_{a}^{\prime },\,{A}_{ab,\boldsymbol{q}}$ and ${A}_{ab,\boldsymbol{q}}^{\prime }$ are respectively detailed in equations (A6b)–(A6e) in the appendix A. The mean-field coefficients ${\mathcal{C}}_{\chi \chi },\,{\mathcal{C}}_{\chi \chi }^{\prime },\,{\mathcal{D}}_{\alpha \beta }$ and ${\mathcal{D}}_{\alpha \beta }^{\prime }$ are functions of temperature and magnon energy through the Bose–Einstein distribution n_{χ, 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, equations (18) and (19) 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 $({\mathcal{K}}_{x}=0)$ and biaxial $({\mathcal{K}}_{x}\ne 0)$ AF systems with hexagonal and square lattice structures that are two common lattice geometries.

In this section, we investigate the effect of nonlinear magnon interactions on magnon dispersion, magnon bandgaps, and AF resonance modes of easy-axis and easy-plane antiferromagnetic insulators.

4.1. Uniaxial easy-axis 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 equation (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 [31, 59], hydrodynamic regime [59] and bulk spin transport [55, 62].

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 equation (11) 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 A1; is given by,

$$\begin{equation}{H}_{\text{MF}}^{\mathrm{u}(4)}=\sum\limits _{\boldsymbol{q}}\left[{\mathcal{C}}_{\boldsymbol{q}}{\alpha }_{\boldsymbol{q}}^{\mathrm{u}{\dagger}}{\alpha }_{\boldsymbol{q}}^{\mathrm{u}}+{\mathcal{C}}_{\boldsymbol{q}}^{\prime }{\beta }_{\boldsymbol{q}}^{\mathrm{u}{\dagger}}{\beta }_{\boldsymbol{q}}^{\mathrm{u}}+{\mathcal{D}}_{\boldsymbol{q}}^{\prime }{\alpha }_{\boldsymbol{q}}^{\mathrm{u}}{\beta }_{-\boldsymbol{q}}^{\mathrm{u}}+\text{H.c.}\right],\end{equation} \tag{ 20 }$$

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

The magnon nonconserving term with strength ${\mathcal{D}}_{\boldsymbol{q}}^{\prime }$, 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,

$$\begin{equation}{\mathcal{C}}_{\boldsymbol{q}}=\frac{1}{2}[{A}_{1}(\vert {u}_{\boldsymbol{q}}^{\mathrm{u}}{\vert }^{2}+\vert {v}_{\boldsymbol{q}}^{\mathrm{u}}{\vert }^{2})-({A}_{2,\boldsymbol{q}}\,{\bar{u}}_{\boldsymbol{q}}^{\mathrm{u}}{\bar{v}}_{\boldsymbol{q}}^{\mathrm{u}}+{\bar{A}}_{2,\boldsymbol{q}}\,{u}_{\boldsymbol{q}}^{\mathrm{u}}{v}_{\boldsymbol{q}}^{\mathrm{u}})],\end{equation} \tag{ 21a }$$

$$\begin{equation}{\mathcal{C}}_{\boldsymbol{q}}^{\prime }=\frac{1}{2}[{A}_{1}(\vert {u}_{\boldsymbol{q}}^{\mathrm{u}}{\vert }^{2}+\vert {v}_{\boldsymbol{q}}^{\mathrm{u}}{\vert }^{2})-({A}_{2,-\boldsymbol{q}}\,{\bar{u}}_{\boldsymbol{q}}^{\mathrm{u}}{\bar{v}}_{\boldsymbol{q}}^{\mathrm{u}}+{\bar{A}}_{2,-\boldsymbol{q}}\,{u}_{\boldsymbol{q}}^{\mathrm{u}}{v}_{\boldsymbol{q}}^{\mathrm{u}})],\end{equation} \tag{ 21b }$$

$$\begin{equation}{\mathcal{D}}_{\boldsymbol{q}}^{\prime }=-2{A}_{1}{\bar{u}}_{\boldsymbol{q}}^{\mathrm{u}}{v}_{\boldsymbol{q}}^{\mathrm{u}}+[{\bar{A}}_{2,\boldsymbol{q}}{({v}_{\boldsymbol{q}}^{\mathrm{u}})}^{2}+{A}_{2,\boldsymbol{q}}{({\bar{u}}_{\boldsymbol{q}}^{\mathrm{u}})}^{2}],\end{equation} \tag{ 21c }$$

the Bogoliubov parameters of the uniaxial case, ${u}_{\boldsymbol{q}}^{\mathrm{u}}$ and ${v}_{\boldsymbol{q}}^{\mathrm{u}}$, and temperature-dependent mean-field coefficients A₁ and A₂ are presented in appendix A. 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}_{\text{m}}^{\mathrm{u}(2)}+{H}_{\text{MF}}^{\mathrm{u}(4)}$. Due to the presence of interband processes, ${\mathcal{D}}^{\prime }$, α and β modes are coupled and the Hamiltonian H^u is no longer diagonal in the {α, β} basis. To diagonalize H^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 }_{\boldsymbol{q}}^{\mu }{\hat{\mu }}_{\boldsymbol{q}}^{{\dagger}}{\hat{\mu }}_{\boldsymbol{q}}+{\varepsilon }_{\boldsymbol{q}}^{\nu }{\hat{\nu }}_{\boldsymbol{q}}^{{\dagger}}{\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 }_{\boldsymbol{q}}^{\mu }={\varepsilon }_{\boldsymbol{q}}^{\nu }\equiv {\varepsilon }_{\boldsymbol{q}}$, where

$$\begin{equation}{\varepsilon }_{\boldsymbol{q}}=\sqrt{{\left({{\epsilon}}_{\boldsymbol{q}}^{\text{u}}+{\mathcal{C}}_{\boldsymbol{q}}\right)}^{2}-\vert {\mathcal{D}}_{\boldsymbol{q}}^{\prime }{\vert }^{2}},\end{equation} \tag{ 22 }$$

with ${{\epsilon}}_{\boldsymbol{q}}^{\text{u}}=zJS\sqrt{{\left(1+2{\mathcal{K}}_{z}/zJ\right)}^{2}-\vert {\gamma }_{\boldsymbol{q}}{\vert }^{2}}$ denotes the dispersion relation of noninteracting magnons in the uniaxial case, i.e., equation (11) with ${\mathcal{K}}_{x}=0$.

The dispersion in equation (22) 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, equations (21a)–(21c), and the eigenenergy ɛ _q , is done self-consistently and involves the calculation of A₁ and A₂ for various temperatures.

In figure 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}}_{\boldsymbol{q}}^{\text{u}}$ at zero temperature (dashed line at figures 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.

Zoom In Zoom Out Reset image size

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

4.2. Biaxial easy-plane 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 equation (11), and the corresponding eigenstates are elliptically polarized magnons [20]. The total interacting magnon Hamiltonian within the mean-field approximation is $H={H}_{\text{m}}^{(2)}+{H}_{\text{MF}}^{(4)}$. The mean-field Hamiltonian (17) can be partitioned into the intra- and interband contribution as ${H}_{\text{MF}}^{(4)}={H}_{\text{intra}}^{(4)}+{H}_{\text{inter}}^{(4)}$, with

$$\begin{equation}{H}_{\text{intra}}^{(4)}=\sum\limits _{\boldsymbol{q}}\left[{\mathcal{C}}_{\alpha \alpha ,\boldsymbol{q}}{\alpha }_{\boldsymbol{q}}^{{\dagger}}{\alpha }_{\boldsymbol{q}}+{\mathcal{C}}_{\beta \beta ,\boldsymbol{q}}{\beta }_{\boldsymbol{q}}^{{\dagger}}{\beta }_{\boldsymbol{q}}+{\mathcal{C}}_{\alpha \alpha ,\boldsymbol{q}}^{\prime }{\alpha }_{\boldsymbol{q}}{\alpha }_{-\boldsymbol{q}}+{\mathcal{C}}_{\beta \beta ,\boldsymbol{q}}^{\prime }{\beta }_{\boldsymbol{q}}{\beta }_{-\boldsymbol{q}}+\text{H.c.}\right],\end{equation} \tag{ 23a }$$

$$\begin{equation}{H}_{\text{inter}}^{(4)}=\sum\limits _{\boldsymbol{q}}\left[{\mathcal{D}}_{\alpha \beta ,\boldsymbol{q}}{\alpha }_{\boldsymbol{q}}^{{\dagger}}{\beta }_{\boldsymbol{q}}+{\mathcal{D}}_{\alpha \beta ,\boldsymbol{q}}^{\prime }{\alpha }_{\boldsymbol{q}}{\beta }_{-\boldsymbol{q}}+\text{H.c.}\right].\end{equation} \tag{ 23b }$$

The intraband contribution of nonlinear magnon scattering consists of diagonal elements, in the α and β operators representation, with coefficients ${\mathcal{C}}_{\chi \chi }$ and the off-diagonal elements with coefficients ${\mathcal{C}}_{\chi \chi }^{\prime }$. ${\mathcal{C}}_{\chi \chi =\alpha \alpha (\beta \beta )}$ renormalizes noninteracting magnon dispersions, obtained within the linear spin-wave theory ${H}_{\text{m}}^{(2)}$, equation (15), while ${\mathcal{C}}_{\chi \chi =\alpha \alpha (\beta \beta )}^{\prime }$ characterize magnon scatterings between states with opposite momenta, q and − q , inside a particular magnon band. The later is equivalent to the annihilation and excitation of an (a) α(β)-magnon pair with zero net linear momentum and spin angular momentum +2ℏ(−2ℏ). The interband nonlinear magnon scatterings are parameterized by off-diagonal terms, ${\mathcal{D}}_{\alpha \beta }$ and ${\mathcal{D}}_{\alpha \beta }^{\prime }$, corresponding to scattering processes of magnons between α and β bands. The total interacting Hamiltonian, H, is no longer diagonalized in the basis of α and β operators. Thus, we do another canonical transformation and find the diagonal basis in terms of new bosonic operators, say μ and ν. The expressions for mean-field parameters are presented in the appendix A2, see equations (A4a)–(A4d). 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 A2, 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 [64]. The importance of this correction is, however, conditioned to small systems, which is a size limit where our mean-field approach might breakdown.

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

$$\begin{equation}{{\epsilon}}_{\boldsymbol{q}}^{\mu }=\sqrt{{({{\epsilon}}_{\boldsymbol{q}}^{\alpha }+{\mathcal{C}}_{\alpha \alpha ,\boldsymbol{q}})}^{2}-\vert {\mathcal{C}}_{\alpha \alpha ,\boldsymbol{q}}^{\prime }{\vert }^{2}},\end{equation} \tag{ 24a }$$

$$\begin{equation}{{\epsilon}}_{\boldsymbol{q}}^{\nu }=\sqrt{{({{\epsilon}}_{\boldsymbol{q}}^{\beta }+{\mathcal{C}}_{\beta \beta ,\boldsymbol{q}})}^{2}-\vert {\mathcal{C}}_{\beta \beta ,\boldsymbol{q}}^{\prime }{\vert }^{2}}.\end{equation} \tag{ 24b }$$

In figure 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 A3. The noninteracting eigenenergies, ${{\epsilon}}_{\boldsymbol{q}}^{\alpha }$ and ${{\epsilon}}_{\boldsymbol{q}}^{\beta }$, 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}}_{0}^{\mu }$ and ${{\epsilon}}_{0}^{\nu }$, shown in the respective insets of panel (a). These figures show that the energy of one magnon band ${{\epsilon}}_{0}^{\nu }$, decreases with temperature, while another band ${{\epsilon}}_{0}^{\mu }$ has a weak temperature dependence for different strengths of the hard-axis anisotropy.

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

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

In this section, we study the effect of nonlinear thermal magnons on the polarization of magnon eigenmodes based on the calculation in appendix A4. In general, magnon polarization is an essential ingredient in magnonic-based information technology [7, 65, 66]. 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 spin 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}_{i}^{x}\rangle }_{\mathrm{max}}/{\langle {S}_{i}^{y}\rangle }_{\mathrm{max}}=1$, see the inset of figure 4 and right panel of figure 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 two AF sublattices [67–71]. 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, quantum magnon squeezing is reduced, but here we show that thermal squeezing becomes dominant at finite temperature [72]. This thermal squeezing in uniaxial AF systems, arising from an effective attractive magnon interaction and breaking of the local U(1) symmetry, see equation (20) and discussions below it.

The resulting temperature-dependent magnonic modes exhibit anisotropic spatial oscillations. This property is evidenced by the ellipticity, which is determined [20] from the Bogolioubov transformation, (equations (6) and (9)), and diagonalization matrix of the Hamiltonian $H={H}_{\text{m}}^{(2)}+{H}_{\text{MF}}^{(4)}$, that in turn relates the spin components and diagonal basis, $\hat{\mu }$ and $\hat{\nu }$, of operators, see appendix A4 for more details. The result is obtained for the uniaxial (e_uni) and biaxial (e_bi) AF square-lattice systems, and presented at figures 4 and 5 as a function of temperature for different magnetic anisotropy strengths. Similar results is also obtained in hexagonal lattices (data not shown).

In figure 4, we plot the ellipticity e_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 equation (20), and thus is termed thermal squeezed magnons [72], unlike its quantum counterpart defined at T = 0 [67–70]. In the biaxial AF system the ellipticity, e_bi, of μ- and ν-magnonic eigenmode, are plotted at figure 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}_{\text{bi}}^{\mu }$ and ${e}_{\text{bi}}^{\nu }$, at q = 0 are represented (left panels) by blue and red curves, respectively. The magnonic μ(ν)-mode becomes more (less) elliptical, being ${e}_{\text{bi}}^{\mu } > 1$ and ${e}_{\text{bi}}^{\nu }< 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}_{\text{bi}}^{\mu }$ and ${e}_{\text{bi}}^{\nu }$, 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 ν mode has lower energy and this band becomes more populated than the μ mode; therefore the nonlinear effects are stronger in this lower band, see figure 2.

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 [73, 74] also leads to the same effective renormalized Hamiltonian within the Hartree–Fock approximation [42, 49]. 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.

MSh acknowledges Dr Morteza Mohseni for fruitful discussions. MSh and JF acknowledge support from Carl Tryggers Stiftelse, and JF 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 Centers of Excellence funding scheme, Project No. 262633, 'QuSpin'.

All data that support the findings of this study are included within the article (and any supplementary files).

In this appendix, we present the treatment of the magnon–magnon Hamiltonian ${H}_{b}^{(4)}$ within the mean-field approximation. We define the following mean-field parameters [60],

$$\begin{equation}m=\frac{2}{NS}\sum\limits _{\boldsymbol{q}}\langle {a}_{\boldsymbol{q}}^{{\dagger}}{a}_{\boldsymbol{q}}\rangle =\frac{2}{NS}\sum\limits _{\boldsymbol{q}}\langle {b}_{\boldsymbol{q}}^{{\dagger}}{b}_{\boldsymbol{q}}\rangle ,\end{equation} \tag{ A1a }$$

$$\begin{equation}{m}^{\prime }=\frac{2}{NS}\sum\limits _{\boldsymbol{q}}{\gamma }_{-\boldsymbol{q}}\langle {a}_{\boldsymbol{q}}{b}_{\boldsymbol{q}}^{{\dagger}}\rangle ,\end{equation} \tag{ A1b }$$

$$\begin{equation}p=\frac{2}{NS}\sum\limits _{\boldsymbol{q}}\langle {a}_{\boldsymbol{q}}{a}_{-\boldsymbol{q}}\rangle =\frac{2}{NS}\sum\limits _{\boldsymbol{q}}\langle {b}_{\boldsymbol{q}}{b}_{-\boldsymbol{q}}\rangle ,\end{equation} \tag{ A1c }$$

$$\begin{equation}{p}^{\prime }=\frac{2}{NS}\sum\limits _{\boldsymbol{q}}{\gamma }_{-\boldsymbol{q}}\langle {a}_{\boldsymbol{q}}{b}_{-\boldsymbol{q}}\rangle ,\end{equation} \tag{ A1d }$$

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

A1. Uniaxial easy-axis AF system

In the uniaxial AF case, where the magnon basis reads as equation (6), the mean-field parameters become,

$$\begin{equation}m=\frac{2}{NS}\sum\limits _{\boldsymbol{q}}\left(\vert {u}_{\boldsymbol{q}}^{u}{\vert }^{2}+\vert {v}_{\boldsymbol{q}}^{u}{\vert }^{2}\right)\,{n}_{\boldsymbol{q}}+\vert {v}_{\boldsymbol{q}}^{u}{\vert }^{2},\end{equation} \tag{ A2a }$$

$$\begin{equation}{m}^{\prime }=p=0,\end{equation} \tag{ A2b }$$

$$\begin{equation}{p}^{\prime }=-\frac{2}{NS}\sum\limits _{\boldsymbol{q}}{\gamma }_{-\boldsymbol{q}}{u}_{\boldsymbol{q}}^{u}{v}_{\boldsymbol{q}}^{u}\,(2{n}_{\boldsymbol{q}}+1).\end{equation} \tag{ A2c }$$

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

$$\begin{equation}{\mathcal{C}}_{\boldsymbol{q}}=\frac{1}{2}[{A}_{1}(\vert {u}_{\boldsymbol{q}}^{u}{\vert }^{2}+\vert {v}_{\boldsymbol{q}}^{u}{\vert }^{2})-({A}_{2,\boldsymbol{q}}\,{\bar{u}}_{\boldsymbol{q}}^{u}{\bar{v}}_{\boldsymbol{q}}^{u}+{\bar{A}}_{2,\boldsymbol{q}}\,{u}_{\boldsymbol{q}}^{u}{v}_{\boldsymbol{q}}^{u})],\end{equation} \tag{ A3a }$$

$$\begin{equation}{\mathcal{C}}_{\boldsymbol{q}}^{\prime }=\frac{1}{2}[{A}_{1}(\vert {u}_{\boldsymbol{q}}^{u}{\vert }^{2}+\vert {v}_{\boldsymbol{q}}^{u}{\vert }^{2})-({A}_{2,-\boldsymbol{q}}\,{\bar{u}}_{\boldsymbol{q}}^{u}{\bar{v}}_{\boldsymbol{q}}^{u}+{\bar{A}}_{2,-\boldsymbol{q}}\,{u}_{\boldsymbol{q}}^{u}{v}_{\boldsymbol{q}}^{u})],\end{equation} \tag{ A3b }$$

$$\begin{equation}{\mathcal{D}}_{\boldsymbol{q}}^{\mathrm{\prime }}=-2{A}_{1}{\bar{u}}_{\boldsymbol{q}}^{u}{v}_{\boldsymbol{q}}^{u}+({\bar{A}}_{2,\boldsymbol{q}}{({v}_{\boldsymbol{q}}^{u})}^{2}+{A}_{2,\boldsymbol{q}}{({\bar{u}}_{\boldsymbol{q}}^{u})}^{2}),\end{equation} \tag{ A3c }$$

where ${A}_{1}=-\left(zJ+4{\mathcal{K}}_{z}\right)Sm-zJS({p}^{\prime }+{\bar{p}}^{\prime })/2$ and ${A}_{2}=-zJS({\bar{p}}^{\prime }+m)\,{\gamma }_{-\boldsymbol{q}}$, and n _q = (n_{α, q} + n_{β, q} )/2. In the uniaxial case and in the absence of magnetic fields and DMIs n_{α, q} = n_{β, q} .

A2. Biaxial easy-plane AF system

In this case ${\mathcal{K}}_{x}\ne 0$, where the magnon basis reads as equation (9), the mean-field parameters become,

$$\begin{equation}m=\frac{2}{NS}\sum\limits _{\boldsymbol{q}}\left[\left(\vert {u}_{\boldsymbol{q}}{\vert }^{2}+\vert {x}_{\boldsymbol{q}}{\vert }^{2}\right)\,{n}_{\alpha ,\boldsymbol{q}}+\left(\vert {v}_{\boldsymbol{q}}{\vert }^{2}+\vert {w}_{\boldsymbol{q}}{\vert }^{2}\right)\,{n}_{\beta ,\boldsymbol{q}}+\vert {x}_{\boldsymbol{q}}{\vert }^{2}+\vert {w}_{\boldsymbol{q}}{\vert }^{2}\right],\end{equation} \tag{ A4a }$$

$$\begin{equation}p=\frac{2}{NS}\sum\limits _{\boldsymbol{q}}\left[{u}_{\boldsymbol{q}}{x}_{\boldsymbol{q}}\,\left(2{n}_{\alpha ,\boldsymbol{q}}+1\right)+{v}_{\boldsymbol{q}}{w}_{\boldsymbol{q}}\,\left(2{n}_{\beta ,\boldsymbol{q}}+1\right)\right],\end{equation} \tag{ A4b }$$

$$\begin{equation}{m}^{\prime }=\frac{2}{NS}\sum\limits _{\boldsymbol{q}}{\gamma }_{-\boldsymbol{q}}\left[-\left(\vert {u}_{\boldsymbol{q}}{\vert }^{2}+\vert {x}_{\boldsymbol{q}}{\vert }^{2}\right)\,{n}_{\alpha ,\boldsymbol{q}}+\left(\vert {v}_{\boldsymbol{q}}{\vert }^{2}+\vert {w}_{\boldsymbol{q}}{\vert }^{2}\right)\,{n}_{\beta ,\boldsymbol{q}}-\vert {u}_{\boldsymbol{q}}{\vert }^{2}+\vert {v}_{\boldsymbol{q}}{\vert }^{2}\right],\end{equation} \tag{ A4c }$$

$$\begin{equation}{p}^{\prime }=\frac{2}{NS}\sum\limits _{\boldsymbol{q}}{\gamma }_{-\boldsymbol{q}}\left[-{u}_{\boldsymbol{q}}{x}_{\boldsymbol{q}}\,\left(2{n}_{\alpha ,\boldsymbol{q}}+1\right)+{v}_{\boldsymbol{q}}{w}_{\boldsymbol{q}}\,\left(2{n}_{\beta ,\boldsymbol{q}}+1\right)\right],\end{equation} \tag{ A4d }$$

which depend on the Bose–Einstein distribution function of magnons at each mode, n_{α, q} and n_{β, 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_{χ, q} → 0, and thus m = 2∑ _q (|x _q |² + |w _q |²)/NS, p = 2∑ _q (u _q x _q + v _q w _q )/NS, m' = 2∑ _q γ_{− q} (−|u _q |² + |v _q |²)/NS and p' = 2∑ _q γ_{− q} (−u _q x _q + v _q w _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 (17) obey,

$$\begin{align}\hfill {\mathcal{C}}_{\alpha \alpha ,\boldsymbol{q}}& =\frac{1}{2}\left[2{A}_{a}(\vert {u}_{\boldsymbol{q}}{\vert }^{2}+\vert {x}_{\boldsymbol{q}}{\vert }^{2})-({A}_{ab,\boldsymbol{q}}^{\prime }+{\bar{A}}_{ab,\boldsymbol{q}}^{\prime })\,\vert {u}_{\boldsymbol{q}}{\vert }^{2}-({A}_{ab,-\boldsymbol{q}}^{\prime }+{\bar{A}}_{ab,-\boldsymbol{q}}^{\prime })\,\vert {x}_{\boldsymbol{q}}{\vert }^{2}\right.\hfill \\ \hfill & \quad \left.+(2{A}_{a}^{\prime }+2{A}_{b}^{\prime }-{A}_{ab,\boldsymbol{q}}-{A}_{ab,-\boldsymbol{q}})\,{u}_{\boldsymbol{q}}{x}_{\boldsymbol{q}}+(2{\bar{A}}_{a}^{\prime }+2{\bar{A}}_{b}^{\prime }-{\bar{A}}_{ab,\boldsymbol{q}}-{\bar{A}}_{ab,-\boldsymbol{q}})\,{\bar{u}}_{\boldsymbol{q}}{\bar{x}}_{\boldsymbol{q}}\right],\hfill \end{align} \tag{ A5a }$$

$$\begin{align}\hfill {\mathcal{C}}_{\beta \beta ,\boldsymbol{q}}& =\frac{1}{2}\left[2{A}_{a}(\vert {v}_{\boldsymbol{q}}{\vert }^{2}+\vert {w}_{\boldsymbol{q}}{\vert }^{2})+({A}_{ab,\boldsymbol{q}}^{\prime }+{\bar{A}}_{ab,\boldsymbol{q}}^{\prime })\,\vert {v}_{\boldsymbol{q}}{\vert }^{2}+({A}_{ab,-\boldsymbol{q}}^{\prime }+{\bar{A}}_{ab,-\boldsymbol{q}}^{\prime })\,\vert {w}_{\boldsymbol{q}}{\vert }^{2}\right.\hfill \\ \hfill & \quad \left.+(2{A}_{a}^{\prime }+2{A}_{b}^{\prime }+{A}_{ab,\boldsymbol{q}}+{A}_{ab,-\boldsymbol{q}})\,{v}_{\boldsymbol{q}}{w}_{\boldsymbol{q}}+(2{\bar{A}}_{a}^{\prime }+2{\bar{A}}_{b}^{\prime }+{\bar{A}}_{ab,\boldsymbol{q}}+{\bar{A}}_{ab,-\boldsymbol{q}})\,{\bar{v}}_{\boldsymbol{q}}{\bar{w}}_{\boldsymbol{q}}\right],\hfill \end{align} \tag{ A5b }$$

$$\begin{equation}{\mathcal{C}}_{\alpha \alpha ,\boldsymbol{q}}^{\prime }=(2{A}_{a}-{A}_{ab,\boldsymbol{q}}^{\prime }-{{\bar{A}}^{\prime }}_{ab,\boldsymbol{q}}){u}_{\boldsymbol{q}}{\bar{x}}_{\boldsymbol{q}}+({A}_{a}^{\prime }+{A}_{b}^{\prime }-{A}_{ab,\boldsymbol{q}}){u}_{\boldsymbol{q}}^{2}+({\bar{A}}_{a}^{\prime }+{\bar{A}}_{b}^{\prime }-{\bar{A}}_{ab,\boldsymbol{q}}){\bar{x}}_{\boldsymbol{q}}^{2},\end{equation} \tag{ A5c }$$

$$\begin{equation}{\mathcal{C}}_{\beta \beta ,\boldsymbol{q}}^{\prime }=(2{A}_{a}+{A}_{ab,\boldsymbol{q}}^{\prime }+{{\bar{A}}^{\prime }}_{ab,\boldsymbol{q}}){v}_{\boldsymbol{q}}{\bar{w}}_{\boldsymbol{q}}+({A}_{a}^{\prime }+{A}_{b}^{\prime }+{A}_{ab,\boldsymbol{q}}){v}_{\boldsymbol{q}}^{2}+({\bar{A}}_{a}^{\prime }+{\bar{A}}_{b}^{\prime }+{\bar{A}}_{ab,\boldsymbol{q}}){\bar{w}}_{\boldsymbol{q}}^{2},\end{equation} \tag{ A5d }$$

$$\begin{align}\hfill {\mathcal{D}}_{\alpha \beta ,\boldsymbol{q}}& =({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}}+\left({A}_{ab,-\boldsymbol{q}}^{\prime }-{{\bar{A}}^{\prime }}_{ab,-\boldsymbol{q}}\right){x}_{\boldsymbol{q}}{\bar{w}}_{\boldsymbol{q}}-\left({A}_{ab,\boldsymbol{q}}^{\prime }-{{\bar{A}}^{\prime }}_{ab,\boldsymbol{q}}\right){v}_{\boldsymbol{q}}{\bar{u}}_{\boldsymbol{q}}\hfill \\ \hfill & \quad +2({A}_{a}^{\prime }-{A}_{b}^{\prime }){v}_{\boldsymbol{q}}{x}_{\boldsymbol{q}}+2({\bar{A}}_{a}^{\prime }-{\bar{A}}_{b}^{\prime }){\bar{u}}_{\boldsymbol{q}}{\bar{w}}_{\boldsymbol{q}},\hfill \end{align} \tag{ A5e }$$

$$\begin{align}\hfill {\mathcal{D}}_{\alpha \beta ,\boldsymbol{q}}^{\prime }& =({\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}}+\left({\bar{A}}_{ab,-\boldsymbol{q}}^{\prime }-{A}_{ab,-\boldsymbol{q}}^{\prime }\right){\bar{x}}_{\boldsymbol{q}}{v}_{\boldsymbol{q}}+\left({A}_{ab,\boldsymbol{q}}^{\prime }-{{\bar{A}}^{\prime }}_{ab,\boldsymbol{q}}\right){u}_{\boldsymbol{q}}{\bar{w}}_{\boldsymbol{q}}\hfill \\ \hfill & \quad +2({A}_{a}^{\prime }-{A}_{b}^{\prime }){u}_{\boldsymbol{q}}{v}_{\boldsymbol{q}}+2({\bar{A}}_{a}^{\prime }-{\bar{A}}_{b}^{\prime }){\bar{x}}_{\boldsymbol{q}}{\bar{w}}_{\boldsymbol{q}},\hfill \end{align} \tag{ A5f }$$

where the following definitions have been introduced

$$\begin{equation}{A}_{a}=-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]\end{equation} \tag{ A6a }$$

$$\begin{equation}{A}_{a}^{\prime }=-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],\end{equation} \tag{ A6b }$$

$$\begin{equation}{A}_{b}^{\prime }=-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],\end{equation} \tag{ A6c }$$

$$\begin{equation}{A}_{ab,\boldsymbol{q}}=-zJS\,{\gamma }_{-\boldsymbol{q}}\,\left({\bar{p}}^{\prime }+m\right)\end{equation} \tag{ A6d }$$

$$\begin{equation}{A}_{ab,\boldsymbol{q}}^{\prime }=-zJS\,{\gamma }_{-\boldsymbol{q}}\,\left({\bar{m}}^{\prime }+\frac{p+\bar{p}}{4}\right).\end{equation} \tag{ A6e }$$

In the cases evaluated in the main text, i.e., in the absence of magnetic field and DMIs, the Bogoliubov coefficients u _q , v _q , w _q and x _q can be defined as real functions, thus the previous relations for the mean-field coefficients are simplified in the following form,

$$\begin{equation}{\mathcal{C}}_{\alpha \alpha ,\boldsymbol{q}}={{\Delta}}_{\boldsymbol{q}}^{-}\vert {u}_{\boldsymbol{q}}{\vert }^{2}+{{\Delta}}_{-\boldsymbol{q}}^{-}\vert {x}_{\boldsymbol{q}}{\vert }^{2}+2\,\text{Re}\left[{u}_{\boldsymbol{q}}{x}_{\boldsymbol{q}}\,\left(\xi -\text{Re}\,{{\Gamma}}_{q}\right)\right],\end{equation} \tag{ A7a }$$

$$\begin{equation}{\mathcal{C}}_{\beta \beta ,\boldsymbol{q}}={{\Delta}}_{\boldsymbol{q}}^{+}\vert {v}_{\boldsymbol{q}}{\vert }^{2}+{{\Delta}}_{\boldsymbol{q}}^{+}\vert {w}_{\boldsymbol{q}}{\vert }^{2}+2\,\text{Re}\left[{v}_{\boldsymbol{q}}{w}_{\boldsymbol{q}}\,\left(\xi +\text{Re}\,{{\Gamma}}_{q}\right)\right],\end{equation} \tag{ A7b }$$

$$\begin{equation}{\mathcal{C}}_{\alpha \alpha ,\boldsymbol{q}}^{\prime }=2{{\Delta}}_{\boldsymbol{q}}^{-}{u}_{\boldsymbol{q}}{x}_{\boldsymbol{q}}+{{\Lambda}}_{\boldsymbol{q}}^{-}{u}_{\boldsymbol{q}}^{2}+{\bar{{\Lambda}}}_{\boldsymbol{q}}^{-}{x}_{\boldsymbol{q}}^{2},\end{equation} \tag{ A7c }$$

$$\begin{equation}{\mathcal{C}}_{\beta \beta ,\boldsymbol{q}}^{\prime }=2{{\Delta}}_{\boldsymbol{q}}^{+}{v}_{\boldsymbol{q}}{w}_{\boldsymbol{q}}+{{\Lambda}}_{\boldsymbol{q}}^{+}{v}_{\boldsymbol{q}}^{2}+{\bar{{\Lambda}}}_{\boldsymbol{q}}^{+}{w}_{\boldsymbol{q}}^{2},\end{equation} \tag{ A7d }$$

$$\begin{equation}{\mathcal{D}}_{\alpha \beta ,\boldsymbol{q}}={{\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}},\end{equation} \tag{ A7e }$$

$$\begin{equation}{\mathcal{D}}_{\alpha \beta ,\boldsymbol{q}}^{\mathrm{\prime }}={\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 }^{\mathrm{\prime }}{u}_{\boldsymbol{q}}{v}_{\boldsymbol{q}}+2\;{\bar{\xi }}^{\mathrm{\prime }}\;{x}_{\boldsymbol{q}}{w}_{\boldsymbol{q}},\end{equation} \tag{ A7f }$$

where we define,

$$\begin{equation}{{\Delta}}_{\boldsymbol{q}}^{\pm }={A}_{a}\pm \text{Re}\left[{A}_{ab,\boldsymbol{q}}^{\prime }\right],\end{equation} \tag{ A8a }$$

$$\begin{equation}\xi ={A}_{a}^{\prime }+{A}_{b}^{\prime },\quad {\xi }^{\prime }={A}_{a}^{\prime }-{A}_{b}^{\prime },\end{equation} \tag{ A8b }$$

$$\begin{equation}{{\Gamma}}_{\boldsymbol{q}}=({A}_{ab,-\boldsymbol{q}}-{A}_{ab,\boldsymbol{q}}),\quad {{\Gamma}}_{\boldsymbol{q}}^{\prime }=({A}_{ab,-\boldsymbol{q}}+{A}_{ab,\boldsymbol{q}}),\end{equation} \tag{ A8c }$$

$$\begin{equation}{{\Lambda}}_{\boldsymbol{q}}^{\pm }=\xi \pm {A}_{ab,\boldsymbol{q}},\end{equation} \tag{ A8d }$$

$$\begin{equation}{{\Theta}}_{\boldsymbol{q}}=({A}_{ab,\boldsymbol{q}}^{\prime }-{\bar{A}}_{ab,\boldsymbol{q}}^{\prime }).\end{equation} \tag{ A8e }$$

A3. Eigenenergies of a biaxial easy-plane AF system

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

$$\begin{align}\hfill 2{\tilde{{\epsilon}}}_{\mu ,\nu }^{2}& =\mathcal{X}\pm \left[{\mathcal{X}}^{2}-\left({\mathcal{D}}_{\alpha \beta }^{\prime 2}\left(0.25{\bar{\mathcal{D}}}_{\alpha \beta }^{\prime 2}-{\bar{\mathcal{C}}}_{\alpha \alpha }^{\prime }{\bar{\mathcal{C}}}_{\beta \beta }^{\prime }\right)+4{\left({{\epsilon}}_{\alpha }+0.5{\mathcal{C}}_{\alpha \alpha }\right)}^{2}\left({\left({{\epsilon}}_{\beta }+0.5{\mathcal{C}}_{\beta \beta }\right)}^{2}-{\mathcal{C}}_{\beta \beta }^{\prime }{\bar{\mathcal{C}}}_{\beta \beta }^{\prime }\right)\right)\right.\hfill \\ \hfill & \quad \left.+{\mathcal{C}}_{\alpha \alpha }^{\prime }\left(4{\left({{\epsilon}}_{\beta }+0.5{\mathcal{C}}_{\beta \beta }\right)}^{2}{\bar{\mathcal{C}}}_{\alpha \alpha }^{\prime }+{\bar{\mathcal{D}}}_{\alpha \beta }^{\prime 2}{\mathcal{C}}_{\beta \beta }^{\prime }-4{\bar{\mathcal{C}}}_{\alpha \alpha }^{\prime }\vert {\mathcal{C}}_{\beta \beta }^{\prime }{\vert }^{2}-2\left({{\epsilon}}_{\beta }+0.5{\mathcal{C}}_{\beta \beta }\right){\bar{\mathcal{D}}}_{\alpha \beta }^{\prime }{\mathcal{D}}_{\alpha \beta }+{\bar{\mathcal{C}}}_{\beta \beta }^{\prime }{\mathcal{D}}_{\alpha \beta }^{2}\right)\right.\hfill \\ \hfill & \quad \left.-2{\mathcal{D}}_{\alpha \beta }^{\prime }{\bar{\mathcal{D}}}_{\alpha \beta }\left(\left({{\epsilon}}_{\beta }+0.5{\mathcal{C}}_{\beta \beta }\right){\bar{\mathcal{C}}}_{\alpha \alpha }^{\prime }-0.25{\bar{\mathcal{D}}}_{\alpha \beta }^{\prime }{\mathcal{D}}_{\alpha \beta }\right)-\left(-{\bar{\mathcal{C}}}_{\alpha \alpha }^{\prime }{\mathcal{C}}_{\beta \beta }^{\prime }+0.25{\mathcal{D}}_{\alpha \beta }^{2}\right){\bar{\mathcal{D}}}_{\alpha \beta }^{2}\right.\hfill \\ \hfill & \quad -2\left({{\epsilon}}_{\alpha }+0.5{\mathcal{C}}_{\alpha \alpha }\right)\left({\mathcal{D}}_{\alpha \beta }^{\prime }{\bar{\mathcal{C}}}_{\beta \beta }^{\prime }{\mathcal{D}}_{\alpha \beta }+{\bar{\mathcal{D}}}_{\alpha \beta }^{\prime }{\mathcal{C}}_{\beta \beta }^{\prime }{\bar{\mathcal{D}}}_{\alpha \beta }-\left({{\epsilon}}_{\beta }+0.5{\mathcal{C}}_{\beta \beta }\right)\left(\vert {\mathcal{D}}_{\alpha \beta }^{\prime }{\vert }^{2}+\vert {\mathcal{D}}_{\alpha \beta }{\vert }^{2}\right)\right){\left.\right]}^{\frac{1}{2}},\hfill \end{align} \tag{ A9 }$$

with $\mathcal{X}={\left({{\epsilon}}_{\alpha }+0.5{\mathcal{C}}_{\alpha \alpha }\right)}^{2}+{\left({{\epsilon}}_{\beta }+0.5{\mathcal{C}}_{\beta \beta }\right)}^{2}+0.5\left({\left\vert {\mathcal{D}}_{\alpha \beta }\right\vert }^{2}-{\left\vert {\mathcal{D}}_{\alpha \beta }^{\prime }\right\vert }^{2}\right)-{\left\vert {\mathcal{C}}_{\alpha \alpha }^{\prime }\right\vert }^{2}-{\left\vert {\mathcal{C}}_{\beta \beta }^{\prime }\right\vert }^{2}$. In square lattices, as we have already discussed in the main text, the interband scattering processes are absent ${\mathcal{D}}_{\alpha \beta }={\mathcal{D}}_{\alpha \beta }^{\prime }=0$, and thus equation (A9) reduces to equation (24) in the main text. The result of equation (A9), is computed for cubic lattices and presented at figure 6 for different temperatures.

Zoom In Zoom Out Reset image size

A4. Ellipticity

The ellipticity is defined as $e={\langle {S}_{i}^{x}\rangle }_{\mathrm{max}}/{\langle {S}_{i}^{y}\rangle }_{\mathrm{max}}$. In the uniaxial easy-axis AF case, the elipticity of each AF eigenmode is given by,

$$\begin{equation}{e}_{\text{uni}}^{\mu }=\frac{\sqrt{A+{\varepsilon }_{\boldsymbol{q}}}+\sqrt{A-{\varepsilon }_{\boldsymbol{q}}}}{\sqrt{A+{\varepsilon }_{\boldsymbol{q}}}-\sqrt{A-{\varepsilon }_{\boldsymbol{q}}}},\end{equation} \tag{ A10a }$$

$$\begin{equation}{e}_{\text{uni}}^{\nu }=\frac{\sqrt{{A}^{\prime }+{\varepsilon }_{\boldsymbol{q}}}+\sqrt{{A}^{\prime }-{\varepsilon }_{\boldsymbol{q}}}}{\sqrt{{A}^{\prime }+{\varepsilon }_{\boldsymbol{q}}}-\sqrt{{A}^{\prime }-{\varepsilon }_{\boldsymbol{q}}}},\end{equation} \tag{ A10b }$$

where $A={{\epsilon}}_{\boldsymbol{q}}^{\text{u}}+0.5{\mathcal{C}}_{\boldsymbol{q}}$ and ${A}^{\prime }={{\epsilon}}_{\boldsymbol{q}}^{\text{u}}+0.5{\mathcal{C}}_{\boldsymbol{q}}^{\prime }$. Since ${\mathcal{C}}_{\boldsymbol{q}}={\mathcal{C}}_{\boldsymbol{q}}^{\prime }$, two modes have the same ellipticity.

In the biaxial easy-plane AF system, we find the following relations for ellipticity of each eigenmode,

$$\begin{equation}{e}_{\text{bi}}^{\mu }=\frac{{P}_{11}+{P}_{13}}{{P}_{11}-{P}_{13}},\end{equation} \tag{ A11a }$$

$$\begin{equation}{e}_{\text{bi}}^{\nu }=\frac{{P}_{14}+{P}_{12}}{{P}_{14}-{P}_{12}},\end{equation} \tag{ A11b }$$

where we have defined the following parameters,

$$\begin{equation}{P}_{11}={u}_{\boldsymbol{q}}{Q}_{11}+{v}_{\boldsymbol{q}}{Q}_{21}+{x}_{\boldsymbol{q}}{Q}_{13}+{w}_{\boldsymbol{q}}{Q}_{23},\end{equation} \tag{ A12a }$$

$$\begin{equation}{P}_{13}={u}_{\boldsymbol{q}}{Q}_{13}+{v}_{\boldsymbol{q}}{Q}_{23}+{x}_{\boldsymbol{q}}{Q}_{11}+{w}_{\boldsymbol{q}}{Q}_{21},\end{equation} \tag{ A12b }$$

$$\begin{equation}{P}_{14}={u}_{\boldsymbol{q}}{Q}_{12}+{v}_{\boldsymbol{q}}{Q}_{22}+{x}_{\boldsymbol{q}}{Q}_{14}+{w}_{\boldsymbol{q}}{Q}_{24},\end{equation} \tag{ A12c }$$

$$\begin{equation}{P}_{12}={u}_{\boldsymbol{q}}{Q}_{14}+{v}_{\boldsymbol{q}}{Q}_{24}+{x}_{\boldsymbol{q}}{Q}_{12}+{w}_{\boldsymbol{q}}{Q}_{22}.\end{equation} \tag{ A12d }$$

The matrix Q is the transformation matrix that diagonalizes the total interacting Hamiltonian $H={H}_{\text{m}}^{(2)}+{H}_{\text{MF}}^{(4)}$, where ${H}_{\text{MF}}^{(4)}={H}_{\text{intra}}^{(4)}+{H}_{\text{inter}}^{(4)}$ is the mean-field interaction, see equations (23a) and (23b) in the main text. In the diagonal basis, the Hamiltonian is written as

$$\begin{equation}H=\frac{1}{2}\sum\limits _{\boldsymbol{q}}{({\Upsilon})}^{{\dagger}}{\Omega}\,({\Upsilon}),\end{equation} \tag{ A13 }$$

where (ϒ) = Q(Ψ) is the new basis in which the Hamiltonian is diagonal $({\Upsilon})={\left({\mu }_{\boldsymbol{q}},{\nu }_{\boldsymbol{q}},{\mu }_{-\boldsymbol{q}}^{{\dagger}},{\nu }_{-\boldsymbol{q}}^{{\dagger}}\right)}^{\mathrm{T}}$, while the old basis is $({\Psi})={\left({\alpha }_{\boldsymbol{q}},{\beta }_{\boldsymbol{q}},{\alpha }_{-\boldsymbol{q}}^{{\dagger}},{\beta }_{-\boldsymbol{q}}^{{\dagger}}\right)}^{\mathrm{T}}$, and Ω is a 4 × 4 diagonal matrix, with elements equivalent to the eigenenergies of the problem and

$$\begin{equation}\mathbf{Q}=\left(\begin{matrix}\hfill {\mathbf{Q}}_{1}\hfill & \hfill {\mathbf{Q}}_{2}\hfill \\ \hfill {\bar{\mathbf{Q}}}_{2}\hfill & \hfill {\bar{\mathbf{Q}}}_{1}\hfill \end{matrix}\right),\qquad {\mathbf{Q}}_{1}=\left(\begin{matrix}\hfill {Q}_{11}\hfill & \hfill {Q}_{12}\hfill \\ \hfill {Q}_{21}\hfill & \hfill {Q}_{22}\hfill \end{matrix}\right),\qquad {\mathbf{Q}}_{2}=\left(\begin{matrix}\hfill {Q}_{13}\hfill & \hfill {Q}_{14}\hfill \\ \hfill {Q}_{23}\hfill & \hfill {Q}_{24}\hfill \end{matrix}\right)\end{equation} \tag{ A14 }$$

being the transformation matrix [20].