The Kuramoto model, describing the synchronization dynamics of coupled oscillators, has been generalized in many ways over the past years. One recent extension of the model replaces the oscillators, originally characterized by a single phase, by particles with D-1 internal phases, represented by a point on the surface of the unit D-sphere. Particles are then more easily represented by D-dimensional unit vectors than by D-1 spherical angles. However, numerical integration of the state equations should ensure that the propagated vectors remain unit and that particles rotate on the sphere as predicted by the dynamical equations. As discussed in [1] integration of the three-dimensional Kuramoto model using Euler's method with time step $Δt$ not only changes the norm of the vectors but produces a small rotation of the particles around the wrong axis. Importantly, the error in the axis' direction does not vanish in the limit $Δt \rightarrow 0$. Therefore, instead of displacing the unit vectors in the direction of the velocity one should performed a sequence of direct small rotations, as dictated by the equations of motion. This keeps the particles on the sphere at all times, ensuring exact norm preservation, and rotates the particles around the proper axis for small $Δt$ [1]. Here I propose an alternative way to do such integration by rotations in 3D that can be generalized to more dimensions using Cayley-Hamilton's theorem. Explicit formulas are provided for 2, 3 and 4 dimensions. I also compare the results with the forth order Runge-Kutta method, which seems to provide accurate results even requiring renormalization of the vectors after each integration step.

18 pages, 3 figures