Runge-Kutta conservation of quadratic invariants ($b_i{a}_{ij}+b_j{a}_{ji}=b_i{b}_j$)

For stage states $Y_i=y_n+h{\sum }_j{a}_{ij}{k}_j$ and a vector field with $Y_i^T{k}_i=0$, expansion of the new squared norm gives $\Vert y_{n+1}{\Vert }^2-\Vert y_n{\Vert }^2=-h^2{\sum }_{i,j}{M}_{ij}{k}_i^T{k}_j$, where $M_{ij}=b_i{a}_{ij}+b_j{a}_{ji}-b_i{b}_j$. Hence $M=0$ preserves the quadratic form exactly for every well-defined step. The same argument with a fixed symmetric matrix in the inner product preserves any quadratic invariant of the vector field. In particular it preserves the Euclidean norm for a skew-symmetric matrix differential equation.