Solution (source code)

= Solution

Use the <Holstein–Primakoff transformation> on the physical <occupation number> subspace $0\le n\le2S$. With the <canonical commutation relation> $[a,a^\dagger]=1$, the ordered square root gives
$$
S^+|n\rangle=\sqrt{n(2S-n+1)}\,|n-1\rangle,\qquad S^-|n\rangle=\sqrt{(n+1)(2S-n)}\,|n+1\rangle.
$$
Both endpoints are respected: $S^+|0\rangle=0$ and $S^-|2S\rangle=0$. Consequently,
$$
[S^+,S^-]|n\rangle=\big[(n+1)(2S-n)-n(2S-n+1)\big]|n\rangle=2(S-n)|n\rangle=2S^z|n\rangle.
$$
Since these <Fock states> form a <basis> of the physical spin space, \b[$[S^+,S^-]=2S^z$] there. Similarly $[S^z,S^\pm]=\pm S^\pm$. The <Holstein–Primakoff occupation constraint> is essential: unrestricted bosonic occupation would not represent a spin-$S$ <Hilbert space>.