First of all, thank you for the quick reply and all the awesome work done in ITensor.
The operator is indeed bosonic so that shouldn't be a problem and I think I understand how to modify the Hubbard site based on what you've said. I will try to implement this and see if it works.
I haven't thought about writing the wavefunction directly as a MPS but it does seem to be a better idea, especially since the entangler hamiltonian is highly non-local. I will take a look at how to do this, thanks for the suggestion.