\( \def\D{\mathcal{D}} \)
\(n\)-dimensional rigid body equations
Angular velocity and kinetic energy
A rigid body rotating with matrix \(R(t)\in O(n)\) has kinetic energy \[T = {\textstyle\frac12}\int\|\dot R x\|^2 dm(x).\] With \(\dot R = R\Omega\) this becomes \[T = {\textstyle\frac12}\mathrm{tr}\bigl(\Omega\D\Omega^T\bigr),\] where \(\D = \int xx^T\;dm(x)\) is the so-called inertia dyadic or inertia tensor, as follows. \begin{eqnarray} \int\|\dot R x\|^2 dm(x) &=& \int(R\Omega x)^T(R\Omega x)\; dm(x)\cr &=& \int x^T\Omega^T\Omega x\; dm(x)\cr &=& \int \mathrm{tr}\bigl(x^T\Omega^T\Omega x\bigr)\; dm(x)\cr &=& \mathrm{tr}\bigl(\Omega \int xx^T\;dm(x)\; \Omega^T \bigr) \cr &=& \mathrm{tr}\bigl(\Omega \D \Omega^T \bigr). \end{eqnarray}
Angular momentum
Instead of \(J=\mathbf{r}\times\mathbf{p}\) we define the angular momentum of a particle moving in \(\mathbb{R}^n\) to be \[J = \mathbf{p}\mathbf{r}^T - \mathbf{r}\mathbf{p}^T\] which is an \(n\times n\) skew-symmetric matrix (here the position vector \(\mathbf{r}\) and the momentum vector \(\mathbf{p}\) are written as column vectors).
This agrees with the usual convention in \(\mathbb{R}^3\), with identification of skew-symmetric matrices and vectors.
In a rigid body with fixed centre of mass, we have \(y=Rx\) and so \(\dot y = \dot R x = R\Omega x\) and the total angular momentum (in space) is \begin{eqnarray} J_s &=& \int (\dot y y^T - y\dot y^T) \; dm(x)\cr &=& \int \bigl(R\Omega xx^T R^T -Rxx^T\Omega^TR^T\bigr) \; dm(x)\cr &=& R \int \bigl(\Omega xx^T -xx^T\Omega^T\bigr) \; dm(x)\;R^T\cr &=& R(\Omega\D - \D\Omega^T)R^T. \end{eqnarray}
The Lagrangian in body coordinates (free rigid body) is \[\mathcal{L} = {\textstyle\frac12}\mathrm{tr}(\Omega\D\Omega^T).\] Defining as usual the conjugate angular momentum by \[\left<\mu,\xi\right> := \left<\D\mathcal{L},\xi\right>\] gives (with the convention \(\left<\mu,\xi\right> = \frac12\mathrm{tr}(\mu^T\xi)\)) \[\mu = \Omega\D-\D\Omega^T.\] Since \(\Omega\) is skew-symmetric, this can be written \[\mu = \Omega\D+\D\Omega.\] It follows that \(J_s = R\mu R^T\) as expected since \(\mu\) is the angular momentum in the body (also denoted \(J_b\)).
Equation of motion
Since \(J_s\) is conserved, differentiating \(J_s=R\mu R^T\) gives \[\dot \mu = \mu\Omega-\Omega\mu = [\mu,\Omega].\]
Principal basis
Choose a basis so that \(\D=\mathrm{diag}[\lambda_1,\dots,\lambda_n]\) is diagonal (where of course \(\lambda_i\geq 0\)). Then, \[\mu_{ij} = \Omega_{ij}(\lambda_i+\lambda_j).\] It follows that the linear map \(\Omega\mapsto\mu\) is invertible iff \(\lambda_i+\lambda_j>0\) (for all \(i\neq j\)). This means that at most one of the \(\lambda_i\) can be zero.
The equation of motion becomes \[\dot\Omega_{ij} = \frac{\lambda_i-\lambda_j}{\lambda_i+\lambda_j}(\Omega^2)_{ij}.\]
Note that if there's some symmetry in the body so that \(\lambda_i=\lambda_j\) then the corresponding \(\Omega_{ij}\) is constant.