Using higher order derivatives of quaternions in equations of motion It is common to look at the orientation of a rigid body in term of a quaternion which encodes an axis and angle with a vector and scalar.
$$ \boldsymbol{q} = \pmatrix{ \boldsymbol{\hat{z}} \sin\left( \tfrac{\theta}{2} \right) \\ \cos \left( \tfrac{\theta}{2} \right)} \tag{1}$$
and then it is common to convert the rotational velocity vector $\boldsymbol{\vec{\omega}}$ into a time derivative of the orientation for use in simulation integrations.
$$ \dot{\boldsymbol{q}} = \tfrac{1}{2} \pmatrix{ \boldsymbol{\vec{\omega}} \\ 0 } \boldsymbol{q} \tag{2} $$
$$ \ddot{\boldsymbol{q}} = \tfrac{1}{2} \pmatrix{ \boldsymbol{\vec{\alpha}} \\ 0} \boldsymbol{q} + \dot{\boldsymbol{q}} \boldsymbol{q}^{-1} \dot{\boldsymbol{q}} \tag{3} $$
What I am asking about is keeping the angular quantities in quaternion form ($\boldsymbol{q}$, $\dot{\boldsymbol{q}}$ and $\ddot{\boldsymbol{q}}$) for direct use in equations of motion and kinematics, by deriving the regular rotational vectors from (2) and its derivative (3)
$$ \pmatrix{\boldsymbol{\vec{\omega}} \\ 0} = 2 \dot{\boldsymbol{q}} \boldsymbol{q}^{-1} \tag{4} $$
$$ \pmatrix{\boldsymbol{\vec{\alpha}} \\ 0} = 2 \left( \ddot{\boldsymbol{q}} - \dot{\boldsymbol{q}} \boldsymbol{q}^{-1} \dot{\boldsymbol{q}} \right) \boldsymbol{q}^{-1} \tag{5} $$
So that angular momentum would be on the form of $$\pmatrix{\boldsymbol{\vec{L}} \\ 0} = \mathcal{I}\, \dot{\boldsymbol{q}} \tag{6}$$
with $\mathcal{I}$ being the appropriate 4×4 form making (6) a matrix vector product type of calculation.
Similarly the rotational equations of motion would be
$$ \pmatrix{ \boldsymbol{\vec{\tau}} \\ 0} = \mathcal{I} \ddot{\boldsymbol{q}} + \underbrace{ \boldsymbol{u} }_{\dot{\boldsymbol{q}}\text{ related terms}} \tag{7} $$
In the spirit of Hamilton, I think there might be some significance to the structure of $\mathcal{I}$, and I was wondering if anyone knows of any prior work related to this line of thinking.
Practically I know the above would simplify the rigid body integrators (by keeping all rotational terms in 4-vectors and allowing linear algebra to do its thing). But I think, this might provide some insight into the behavior of rigid body orientation, just as the profile of translational accelerations gives us insight into velocities and displacements.
 A: If you start with the quaternion vector
$$\vec{z}=\begin{bmatrix}
  a(t) \\
  b(t) \\
  c(t) \\
  d(t) \\
\end{bmatrix}=
\begin{bmatrix}
  a \\
  \vec{w} \\
\end{bmatrix}$$
where $\vec{z}\cdot\vec{z}=1$
you obtain:
$$\begin{bmatrix}
  0 \\
  \vec{\omega} \\
\end{bmatrix}=2\,  \underbrace{\begin{bmatrix}
   a & \vec{w}^T \\
   -\vec{w} & a\,I_3+\tilde{w} \\
 \end{bmatrix}}_{Q (4\times 4)}\,\underbrace{\begin{bmatrix}
  \dot{a} \\
  \vec{\dot w} \\
\end{bmatrix}}_{\vec{\dot{z}}}\tag 1$$
where  $Q^T\,Q=I_4$, thus $Q$  matrix is orthogonal matrix , $\vec{\omega}$ is the angular velocity vector and $I_3$ is a $3\times 3$ unity matrix.
$$\tilde{w}=\left[ \begin {array}{ccc} 0&-w_{{z}}&w_{{y}}\\  w_{
{z}}&0&-w_{{x}}\\  -w_{{y}}&w_{{x}}&0\end {array}
 \right] 
$$.
for numerical simulation you can obtain from equation (1)
$$\vec{\dot{z}}=\frac 12 \,\begin{bmatrix}
    0 & -\vec{\omega}^T \\
    \vec{\omega} & -\tilde{\omega} \\
  \end{bmatrix}\,\vec{z}+\frac{P}{2}\,(1-\vec{z}\,\cdot{\vec{z}})\,\vec{z}\tag 2$$
where the second part of the RHS  is to fulfill the requirement that $\vec{z}\cdot\vec{z}$ must be equal to one.
Equation (2) together with the Euler equations are the EOM's, the separation between the kinematic equation (2) and the Euler equation, make the EOM's simple
Of course you can use equation (1) in the Euler equations, but I don't see how you get the requirement that $\vec{z}\cdot\vec{z}$ must be equal to one ?.
The Euler Equations:
$I\,\vec{\dot \omega}+\vec{\omega}\,\times (I\,\vec{ \omega})=\vec{\tau}$
