As discussed in 2020@Huber and Welshman@2013 , the time derivative of a rotation matrix $R$ with an angular velocity $\omega$ is given by $$ \dot{R} = [\omega]_{\times} R \label{eq:r-dot} $$
When rotating with a constant angular velocity, i.e., $\omega \neq 0, \alpha = 0$, for a period of time $\Delta_T$, the new rotation matrix, i.e., new orientation, can be approximated by $$ R(t + \Delta_T) \approx e^{[\hat{\omega}]_{\times} \theta} R(t) \label{eq:r-next} $$ where $[\cdot]_{\times}$ denotes the skew symmetric operator over $\mathbb{R}^3$ vectors, and $\theta = ||\omega|| \Delta_T$, $\hat{\omega} = \frac{\omega}{||\omega||}$ are the the rotation angle and axis; implementation of such computation is provided in Eigen ยท Rotation Logarithm and Exponential.
However, when angular acceleration is not zero, i.e., $\omega \neq 0, \alpha \neq 0$, approximation \eqref{eq:r-next} is not accurate enough, a better solution was proposed by Wilhelm Magnus which reads $$ R(t + \Delta_T) = e^{[\Omega]_{\times}} R(t) $$ where $\Omega \in \mathbb{R}^3$ is the Magnus Expansion and can be calculated by Magnus Series $$ \Omega = \Omega_1 + \Omega_2 + \Omega_3 + \cdots $$
Assuming constant angular acceleration, i.e., $\omega \neq 0, \alpha \neq 0, j = 0$, the summation of the first three terms guarantees that the result converges as long as the angle rotated per time-step is less than $\pi/\sqrt{2}$ radians, or 127 degrees; and they can be calculated as follow. $$ \begin{align} \Omega_1 &= \omega \Delta_T + \frac{1}{2} \alpha \Delta_T^2 \\ \Omega_2 &= \frac{1}{12}[\alpha]_{\times} \omega \Delta_T^3 \\ \Omega_3 &= \frac{1}{240} [\alpha]_{\times} [\alpha]_{\times} \omega \Delta_T^5 \\ \end{align} $$
Similarly in quaternion space, one gets $$ \mathcal{Q}(t + \Delta_T) \approx \frac{1}{2} [\omega(t)]_{\mathcal{Q}} \, \mathcal{Q}(t) $$ and $$ \mathcal{Q}(t + \Delta_T) = e^{ \frac{1}{2} \big[\Omega\big]_{\mathcal{Q}} } \mathcal{Q}(t) $$ where $\big[\cdot\big]_{\mathcal{Q}}$ denotes the pure quaternion operator over $\mathbb{R}^3$ vectors.