Spin-Spin Hamiltonian in coupled harmonic oscillator I was reading about identical particles and i came across this example:

Consider two electrons with spin 1/2.
The Hamiltonian for this system is:
$$Η=\frac{p_1^2}{2m}+\frac{p_2^2}{2m}+\frac{1}{2}m\omega x_1^2+\frac{1}{2}m\omega x_2^2+g\vec{S_1}\cdot\vec{S_2}(x_1-x_2)^2$$

Im not sure what the last term means exactly, is it correct writing in this form?
$$Η=1\otimes1\otimes\frac{p_1^2}{2m}+1\otimes1\otimes\frac{p_2^2}{2m}+1\otimes1\otimes\frac{1}{2}m\omega x_1^2+1\otimes1\otimes\frac{1}{2}m\omega x_2^2+g\vec{S_1}\otimes\vec{S_2}\otimes(x_1-x_2)^2$$
If this is correct i tried changing the coordinates by using this transformation:
$$R=\frac{x_1+x_2}{m},r=x_1-x_2, \mu=\frac{m}{2}, M=2m$$
and then the new Hamiltonian becomes:
$$Η=1\otimes1\otimes\frac{p_r^2}{2\mu}+1\otimes1\otimes\frac{p_R^2}{2M}+1\otimes1\otimes\frac{1}{2}M\omega R^2+1\otimes1\otimes\frac{1}{2}\mu\omega r^2+g\vec{S_1}\otimes\vec{S_2}\otimes(x_1-x_2)^2$$
Does the spin part transform with the new coordinates i set? If so how? Another problem i have is to find the eigenstates of this hamiltonian, the only way i can think of is to consider:
$$Η_0=1\otimes1\otimes\frac{p_r^2}{2\mu}+1\otimes1\otimes\frac{p_R^2}{2M}+1\otimes1\otimes\frac{1}{2}M\omega R^2+1\otimes1\otimes\frac{1}{2}\mu\omega r^2$$ as my original hamiltonian and
$$H_p=g\vec{S_1}\otimes\vec{S_2}\otimes(x_1-x_2)^2$$ but does this hold for any value of $g$?
Last question is there any way to compute the exact solutions to this hamiltonian?
 A: We can express the Hilbert space of this system as $\mathcal{H}_{tot}= \mathcal{H}_{0}\otimes\mathcal{H}_\frac{1}{2}$, where $\mathcal{H}_0$ is the spinless part of the Hilbert space and $\mathcal{H}_{\frac{1}{2}}$ is the spin part. So the answer to your first question, is yes, we can split it up.
Consider the state $\lvert s_1s_2 \rangle$, where  $s_1$ and $s_2$ are the eigenvalues of $\hat{S}^z_1$ and $\hat{S}^z_2$, respectively. Now, the operator $\vec{S}^1 \cdot \vec{S}^2$ is not diagonal in the $\hat{S}^z$ basis. One can diagonalize the $\vec{S}^1 \cdot \vec{S}^2$ operator independently of the rest of the Hamiltonian. I'll leave that for you to do. The point is that we get the eigenvalues $\left\lbrace \frac{1}{4}, \frac{1}{4},\frac{1}{4}, \frac{-3}{4} \right\rbrace$. Let's denote one of these spin eigenvalues as $s_{12}$. Then we can consider a general Hamiltonian eigenstate $\lvert \{n_i\}  \rangle \lvert s_{12} \rangle$.
Acting with the Hamiltonian on this state, we have
$$H \lvert \{n_i\}\rangle \lvert s_{12} \rangle = \left(\frac{p_1^2}{2m} +\frac{p_2^2}{2m} + \frac{1}{2} m \omega x_1^2 + \frac{1}{2} m \omega x_2^2 + gs_{12} (x_1-x_2)^2\right) \lvert \{n_i\}\rangle \lvert s_{12} \rangle.$$
Now we can perform the coordinate transformation, which you state in your question $$R= \frac{x_1+x_2}{m},\quad r=x_1-x_2,\quad \mu=\frac{m}{2},\quad M=2m.$$
This leads to the reduced Hamiltonian on the spinless Hilbert space
$$\mathcal{H}_{s_1s_2}= \frac{p_r^2}{2\mu} +\frac{p_R^2}{2M}+ \frac{1}{2}M \omega R^2 + \frac{1}{2} \mu \omega \left(1 + \frac{2g s_{12}}{\mu \omega} \right)r^2.$$
This is the Hamiltonian of two decoupled harmonic oscillators. Hence the eigenstates $\lvert n_r n_R s_{12} \rangle$ are classified by 3 quantum numbers and the energies of the harmonic oscillators are
$$E_r = \hbar \omega \left(n_r + \frac{1}{2}\right) \text{ and } E_R =  \hbar \omega \left( 1 + \frac{2 g s_{12}}{\mu \omega}\right) \left(n_R + \frac{1}{2}\right).$$
I hope this also answers your last question about computing the exact solutions to the Hamiltonian.
