Solution (source code)

= Solution

Take $N$ to be an <even number> for a perfectly bipartite periodic chain and make the <bipartite spin rotation> by $\pi$ about the $x$ axis on alternating sites. In the local frame the <Néel state> has all spins up, while a bond becomes
$$
\mathbf S_i\cdot\mathbf S_j=-\widetilde S_i^z\widetilde S_j^z+\tfrac12(\widetilde S_i^+\widetilde S_j^++\widetilde S_i^-\widetilde S_j^-).
$$
The <linear spin-wave approximation> keeps $\widetilde S^z=S-a^\dagger a$ and $\widetilde S^\pm\simeq\sqrt{2S}(a,a^\dagger)$. Therefore
$$
H=-NJS^2+JS\sum_{\langle ij\rangle}(n_i+n_j+a_ia_j+a_i^\dagger a_j^\dagger)+O(S^0).
$$
Every site has two neighbours. <Fourier transform> gives $\sum_{\langle ij\rangle}(n_i+n_j)=2\sum_kn_k$ and pairing coefficient $\gamma_k=\cos k$, so
$$
H=-NJS^2+2JS\sum_ka_k^\dagger a_k+JS\sum_k\gamma_k(a_ka_{-k}+a_k^\dagger a_{-k}^\dagger)+O(S^0).
$$
The <bosonic Nambu normal-ordering shift> follows from $a_{-k}a_{-k}^\dagger=1+a_{-k}^\dagger a_{-k}$. It adds $NJS$ inside the matrix expression, which must be subtracted in its constant. Hence
$$
\boxed{H=-NJS(S+1)+JS\sum_k(a_k^\dagger,a_{-k})
\begin{pmatrix}1&\gamma_k\\\gamma_k&1\end{pmatrix}
\begin{pmatrix}a_k\\a_{-k}^\dagger\end{pmatrix}+O(S^0).}
$$
A periodic chain with an <odd number> of sites is frustrated at the boundary and lacks this exact two-sublattice reference; the bulk thermodynamic calculation uses the even-chain sequence.