Sorry about the slow reply - I agree with Jin's response (Jin: feel free to post your comments as a reply instead of just comments if you are confident about them). So I don't think the issue has to do with complex numbers per se but rather the issue of using a non-IQTensor operator in an IQTensor / quantum number conserving context, which causes the code to attempt to compute the divergence of a mixedIQTensor which is not well defined. (Part of the the reason mixedIQTensors exist, in fact.)
A bit more backgroun on why there's such a thing as a mixed IQTensor: earlier on we decided on a design where there is a single SiteSet type for both ITensors and IQTensors for the same type of Hilbert space. (So like SpinHalf, SpinOne, etc. can be used in both an ITensor and IQTensor context.) Thus the return type of the .op method is IQTensor, since IQTensors can always be converted to ITensors. However, there are some operators that are not well-defined IQTensors, such as Sx. So I made a storage type "mixedIQTensor" that can sort of "smuggle" an ITensor in an IQTensor guise.
In the future we may just have two versions of every site set, so like IQSpinHalf as well as SpinHalf.
Another fix might be to make IQTensors less strict about having a well defined divergence as you suggest in the comments. However there's at least one big issue with this: for efficiency when constructing an IQTensor we use the flux, as soon as it can be determined, to work out all of the blocks that can possibly be non-zero and allocate storage for them in one go. That way we only call the allocator once and don't have to resize the memory or have non-contiguous IQTensor memory. Those things could cause big slowdowns.
So in the end the situation as I see it is: (1) one can always work with operators such as S+ and S- that change quantum numbers in a well-defined way for symmetric Hamiltonians; and (2) restricting IQTensors to just these operator types leads to big efficiencies so it's worth it.