Solution (source code)

= Solution

The <Holstein–Primakoff transformation> acts on the physical number-state subspace $0\leq n\leq2S$. Its operator ordering matters: $S^+=\sqrt{2S-\hat n}\,a$ lowers the number first. Therefore
$$
S^+|n\rangle=\sqrt{n(2S-n+1)}|n-1\rangle,\qquad
S^-|n\rangle=\sqrt{(n+1)(2S-n)}|n+1\rangle,
$$
while $S^z|n\rangle=(S-n)|n\rangle$. 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\rangle$ are $(n+1)(2S-n)$ and $n(2S-n+1)$, respectively, so their difference is $2(S-n)$. Thus
$$
\boxed{[S^+,S^-]=2S^z,\qquad[S^z,S^\pm]=\pm S^\pm.}
$$
Also $(S^z)^2+\frac12(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.