Hamiltonian with position-spin coupling I am solving a Hamiltonian including a term $(x\cdot S)^2$.
The Hamiltonian is like this form:
\begin{equation}
H=L\cdot S+(x\cdot S)^2
\end{equation}
where $x$ is the position operator, $L$ is angular momentum operator, and $S$ is spin operator. The eigenvalue for $L^2$ and $S^2$ are $l(l+1)$ and $s(s+1)$. 
If the Hamiltonian only has the first term, it is just spin orbital coupling and it is easy to solve. The total $J=L+S$, $L^2$ and $S^2$ are quantum number. However, when we consider the second term position and spin coupling $(x\cdot S)^2$, it becomes much harder. The total $J$ is still a quantum number. We have $[(x\cdot S)^2, J]=0$. However, $[(x\cdot S)^2,L^2]≠0$, $L$ is not a quantum number anymore.
Anybody have ideas on how to solve this Hamiltonian?
 A: This problem appears interesting for the following reason. Let us write it down in Cartesian coordnates:
$$-\frac{1}{2}\left(\frac{\partial^2\psi}{\partial x^2}+\frac{\partial^2\psi}{\partial y^2}+\frac{\partial^2\psi}{\partial z^2}\right)+\frac{1}{2}(x\cdot S)^2\psi+L\cdot S\psi=E\psi$$
where I have introduced a 1/2 factor for later convenience. Now, I concentrate on x and I consider the operator
$$-\frac{1}{2}\frac{\partial^2}{\partial x^2}+\frac{1}{2}(x\cdot S)^2$$
One can introduce creation and annihilation operators in a similar way as for the harmonic oscillator
$$A_S=\frac{1}{\sqrt{2}}\left(\frac{\partial}{\partial x}+xS\right)$$
and the corresponding eigenvectors will be labeled as $|n,S\rangle$. The next step is to write down $L\cdot S=\frac{1}{2}(J^2-L^2-S^2)$ and we can restate this problem in the form
$$\left(A_S^\dagger A_S+\frac{1}{2}\right)\psi-\frac{1}{2}\left(\frac{\partial^2\psi}{\partial y^2}+\frac{\partial^2\psi}{\partial z^2}\right)+\frac{1}{2}(J^2-L^2-S^2)\psi=E\psi$$
A: I would solve it using a matrix representation. 
If we multiply Pauli matrices by $\frac{\hbar}{2}$ we can work in the following basis:
$$|n;S_z = +\rangle, |n ; S_z = -\rangle$$
Note that 
$$S\cdot L = S_xL_x +S_yL_y +S_zL_z$$
$$X\cdot S = XS_x$$
$$[L,S] =0$$
You get some matrix in the $S_z$ base ($2\times 2$) and diagonalize it (finding eigenstates and eigenvalue).
A: I would recommend to start from inspecting your Hamiltonian symmetry. It is easy to note that it has rotational symmetry around $x$ axis. Thus, the states may be classified in accordance with $J_x$. 
Probably, you could start from Hamiltonian without $xS$ term, write the solution in terms of $J$ and (!) $J_x$ and then examine $(xS)^2$ operator in this basis. 
Maybe, it is more convenient to write your additional term as $(zS)$ which makes it easy to use standard basis from the textbooks.
