I just tried the normalization method and it works surprisingly well, considering that the solver doesn’t know about the norm constraint… I guess the dynamics constraint is helping a lot with that, and, like you said, better integrators might improve it even further.
Many thanks for sharing your solution! I’m glad my code works properly now at least.