Most general form of a spin rotation invariant Hamiltonian? I am told that the most general form of a spin rotation invariant Hamiltonian for two systems 1 and 2 both with spin $S$, i.e., the spin operators
\begin{align}
(\hat{S}_1^x)^2 +(\hat{S}_1^y)^2 + (\hat{S}_1^z)^2 = (\hat{S}_2^x)^2 +(\hat{S}_2^y)^2 + (\hat{S}_2^z)^2 = S(S+1)\hbar^2
\end{align}
is given by
\begin{equation}
\mathcal{H} = \sum_{j=0}^{2S} a_j \bigg(\frac{\mathbf{\hat{S}_1}\cdot\mathbf{\hat{S}_2}}{\hbar}\bigg)^j
\end{equation}
I understand that it should be a function of $\mathbf{\hat{S}_1}\cdot\mathbf{\hat{S}_2}$ but I do not understand why should the sum terminate at $j=2S$. Can someone explain it to me. Thank you.
 A: $\newcommand{\bm}[1]{\mathbf{#1}}$
You need to look at this in terms of spin quantum numbers (i.e., eigenvalues).
$(\bm S_1+\bm S_2)$ can take values $S_\mathrm{tot} = 0,\dots,2S$. Now if we restrict to the subspace with total spin $S_\mathrm{tot}$, we have that
\begin{equation}
\begin{aligned}
2S_\mathrm{tot}(2S_\mathrm{tot}+1) 
= 
(\bm S_1+\bm S_2)^2 &= \bm S_1\cdot \bm S_1 + \bm S_2\cdot \bm S_2 +
2\, \bm S_1\cdot \bm S_2
\\
&= S(S+1) + S(S+1) + 2\,\bm S_1\cdot \bm S_2\ ,
\end{aligned}
\end{equation}
and thus
$$\bm S_1\cdot \bm S_2 = S_\mathrm{tot}(2S_\mathrm{tot}+1)-S(S+1)
\tag{1}
$$
can take $2S+1$ possible values. (Note that this means that $\bm S_1\cdot \bm S_2$ and $\bm S_1+\bm S_2$ are diagonal in the same basis, that is, we can reason
about them as if they were just numbers which can take the corresponding set of values.)
A $\mathrm{SU}(2)$ invariant Hamiltonian of the two spins will take a different value of each subspace of total spin $S_\mathrm{tot}$, i.e., it is of the form (with $\Pi_{S_\mathrm{tot}}$ the projector onto the subspace with total spin $S_\mathrm{tot}$)
$$
\mathcal H = \sum_{S_\mathrm{tot}=0}^{2S} E_{S_\mathrm{tot}} \Pi_{S_\mathrm{tot}}\ .
\tag{2}
$$
Since there is a one-to-one relation between the total spin and the value of $\bm S_1\cdot \bm S_2$ -- Eq (1) --, each projector $\Pi_{S_\mathrm{tot}}$ can be expressed as a function of $\bm S_1\cdot \bm S_2$:
$$\Pi_{S_\mathrm{tot}}=f_{S_\mathrm{tot}}(\bm S_1\cdot \bm S_2)\ .
\tag{3}
$$
This function must be $f_{S_\mathrm{tot}}(\bm S_1\cdot \bm S_2)=1$ for the desired $S_\mathrm{tot}$ (using (1)), and $f(\bm S_1\cdot \bm S_2)=0$ for all other values the total spin can take (again using (1)). This means that we only need to fix the value of $f$ at $2S+1$ points, and thus, it can be choosen to be a polynomial of degree $2S$:
$$
f_{S_\mathrm{tot}}(\bm S_1\cdot \bm S_2) = \sum_{j=0}^{2S} a_{j,S_\mathrm{tot}} (\bm S_1\cdot \bm S_2)^j\ .
$$
Substituting this into (3) and then into (2) gives that
$$
\mathcal H = \sum_{j=0}^{2S} a_{j} (\bm S_1\cdot \bm S_2)^j\ .
$$
