 Solution

The Holstein–Primakoff transformation acts on the physical number-state subspace $0\le n\le 2S$. Its operator ordering matters: $S^{+}=\sqrt{2S-\hat{n}}a$ lowers the number first. Thereforewhile $S^z|n⟩=(S-n)|n⟩$. The square roots vanish at the appropriate ends, so the spin operators preserve the Holstein–Primakoff occupation constraint. The bosonic Fock space outside this subspace is not part of the spin-$S$ representation.

These matrix elements give $[S^z,S^{\pm }]=\pm S^{\pm }$. The two products on $|n⟩$ are $(n+1)(2S-n)$ and $n(2S-n+1)$, respectively, so their difference is $2(S-n)$. ThusAlso $(S^z{)}^2+\frac{1}{2}(S^{+}{S}^{-}+S^{-}{S}^{+})=S(S+1)I$. The transformation therefore has both the correct algebra and the correct spin Casimir operator, not just its small-occupation approximation.