whenever there is no anisotropy and (ii) $\Delta$ does not depend on the type nor the strength of the intralayer coupling.

For evaluating the magnon specific heat (or other thermodynamical quantities), an explicit expression of $J(k)$ is required. In the absence of  such an expression and for low-temerpature range, it is a common practice to assume a long-wave limit. Here, we restricted the expansion of  $J(k)$ to the nearest neighbors only, leading to: $$\hbar\omega_{k}\approx\sqrt{\Delta^{2}+C_{x}k_{x}^{2}+C_{y}k_{y}^{2}+C_{z}k_{z}^{2}}\text{ } \label{dispersion-quadratic}$$ where $C_{x}$ ($\approx C_{y}$ for weak orthorhombic distortion) and $C_{z}$ are functions of the exchange couplings and geometrical factors ($a$ and $c$ are unit-cell parameters):$$\begin{aligned}
C_{x}  &  =16(J_{\bot}+J_{\Vert})J_{\Vert}S^{2}a^{2}+2J_{\bot}S(g\mu_{B}H_{a})a^{2}\label{C-definition}\\
C_{z}  &  =16\ J_{\Vert}^{2}S^{2}c^{2}\nonumber\end{aligned}$$ Then, the zero-field magnon specific heat is (rewritten so as to conform with the notation of Ref.[@Eisele-Keffer-SpinWave]): $$C_{M}(T)=3^{\frac{3}{2}}R\left(  \Delta/\theta\right)  ^{3}\left(
\Delta/2T\pi^{2}\right)  \sum_{m=1}^{\infty}\left[  BesselK(2,m\Delta
/T)+BesselK(4,m\Delta/T)\right]  \label{Cm-bessel}$$ where $BesselK$ represents the modified Bessel function and $$\theta=z\left\vert J_{eff}\right\vert S=\sqrt[3]{3^{\frac{3}{2}}2(C_{x}C_{y}C_{z})^{\frac{1}{2}}/(a^{2}c)} \label{Theta}$$ is a characteristic temperature, based on which $\left\vert J_{eff}\right\vert
$ can be defined as being an effective exchange interaction that couples the magnetic moment to its $z$ nearest neighbors. Notice that both Eq.\[Cm-bessel\] and Eq.\[Theta\] account as well for the weakly orthrohmic-distorted state.

For $T<\Delta$, Eq.\[Cm-bessel\] reduces to the exponential form: $$C_{M}(T)\simeq2\frac{3^{\frac{3}{2}}R\Delta^{\frac{7}{2}}}{\pi^{\frac{3}{2}}\theta^{3}T^{\frac{1}{2}}}\exp(-\Delta/T) \label{Cm-AF-lowT}$$ while for 