The DM interaction has three coordinate-specific terms when splitting it up. Two of these, the DM-x and DM-z terms, are imaginary when we transform them into series of raising and lowering operators. How do we evaluate this? Do we just ignore the i?
\begin{equation} \mathbf{D}_z \cdot (\mathbf{S}_i \times \mathbf{S}_j) = \frac{iD_z}{2} (S^+_i S^-_j - S^-_i S^+_j) \end{equation}
The term you wrote corresponds to a 4x4 matrix \begin{equation} \begin{pmatrix} 0&0&0&0\\ 0&0&-iD_z/2&0\\ 0&iD_z/2&0&0\\ 0&0&0&0 \end{pmatrix} \end{equation} Then, it should become clear that you cannot just ignore the $i$ because without it the matrix would not be Hermitian. You can play around with this matrix to see what kind of states it favors.
While it may be hard to construct intuition on what this term is physically doing, a hint from the classical analogue is that it should be favoring a particular chirality among the spins.
I am not sure what you exactly mean by evaluating. If you want to obtain the expectation value of this term for any particular state $|\psi\rangle$, then you can just calculate $\langle \psi| D_z\cdot(S_i\times S_j) |\psi\rangle$ with the $i$. It will never give you a complex number for that.
