Given two skew lines, $\mathbf{x}_1 + \mathbf{u}_1 t$ and $\mathbf{x}_2 + \mathbf{u}_2 t$ where $\mathbf{u}_1$ and $\mathbf{u}_2$ are unit vectors along the lines, we might be able to rotate each line with angular velocity $\mathbf{w}_1$ and $\mathbf{w}_2$. I am interested in calculating the change of angle between the two lines due to the rotations. The sine of the angle should be $| \mathbf{u}_1 \times \mathbf{u}_2 |$. For simplicity, I dealt with $(\mathbf{u}_1 \times \mathbf{u}_2)^2 = (\mathbf{u}_1 \times \mathbf{u}_2) \cdot (\mathbf{u}_1 \times \mathbf{u}_2)$. Diferentiating this expression with respect to $t$, I get
$$ \frac{d}{dt}(\mathbf{u}_1 \times \mathbf{u}_2)^2 = 2(\mathbf{u}_1 \times \mathbf{u}_2)\cdot (\frac{d\mathbf{u}_1}{dt}\times \mathbf{u}_2 + \mathbf{u}_1 \times \frac{d\mathbf{u}_2}{dt}). $$
Using $\frac{d \mathbf{u}}{d t} = \mathbf{w} \times \mathbf{u}$, the above becomes $$ \frac{d}{dt}(\mathbf{u}_1 \times \mathbf{u}_2)^2 = 2(\mathbf{u}_1 \times \mathbf{u}_2)\cdot ((\mathbf{w}_1\times \mathbf{u}_1)\times \mathbf{u}_2 + \mathbf{u}_1 \times (\mathbf{w}_2 \times \mathbf{u}_2)). $$
Now, using $(\mathbf{A}\times\mathbf{B})\times\mathbf{C} = -(\mathbf{C}\cdot\mathbf{B})\mathbf{A} + (\mathbf{C}\cdot\mathbf{A})\mathbf{B}$ and $\mathbf{A}\times(\mathbf{B}\times\mathbf{C}) = (\mathbf{A}\cdot\mathbf{C})\mathbf{B} - (\mathbf{A}\cdot\mathbf{B})\mathbf{C}$, the vector triple product terms become
$$ (\mathbf{w}_1 \times \mathbf{u}_1) \times \mathbf{u}_2 = -(\mathbf{u}_2\cdot \mathbf{u}_1)\mathbf{w}_1 + (\mathbf{u}_2 \cdot \mathbf{w}_1) \mathbf{u}_1 \\ \mathbf{u}_1 \times (\mathbf{w}_2 \times \mathbf{u}_2) = (\mathbf{u}_1\cdot \mathbf{u}_2)\mathbf{w}_2 - (\mathbf{u}_1 \cdot \mathbf{w}_2) \mathbf{u}_2 $$
then, the above differentiation becomes quite a simple expression:
$$ \frac{d}{dt}(\mathbf{u}_1 \times \mathbf{u}_2)^2 = 2(\mathbf{u}_1 \cdot \mathbf{u}_2)(\mathbf{w}_2 - \mathbf{w}_1)\cdot(\mathbf{u}_1 \times \mathbf{u}_2). $$
because $\mathbf{u}_i \cdot (\mathbf{u}_1 \times \mathbf{u}_2) = 0$ with $i = 1,2$.
The interpretation of this result is tricky for me because of $\mathbf{u}_1 \cdot \mathbf{u}_2$ term. Does it really mean that if $\mathbf{u}_1$ and $\mathbf{u}_2$ were orthogonal initially, the derivative of its cross product is zero regardless of how each line rotates. Or were there error in my calculation so far?
You derivation is correct but the function you're differentiating is the square of the sine of the angle, i.e.
$f(t) = \sin^2 \theta(t) $
So when you differentiate with respect to $t$, you will get
$ f'(t) = 2 \sin \theta \cos \theta \dfrac{d \theta}{d t } $
And this will be $0$ if $\theta = 90^\circ$ because of the $\cos \theta$ term.
An alternative to your derivation is to start with
$ \cos \theta = u_1 \cdot u_2 $
Then
$ (- \sin \theta ) \dfrac{d \theta} {dt} = u_1' \cdot u_2 + u_2' \cdot u_1 $
Since $u_1' = \omega_1 \times u_1 $ and $u_2' = \omega_2 \times u_2 $ , then
$ (- \sin \theta ) \dfrac{d \theta} {dt} = (\omega_1 \times u_1) \cdot u_2 + (\omega_2 \times u_2) \cdot u_1 $
Now, using $ a \cdot (b \times c) = c \cdot ( a \times b) = b \cdot (c \times a) $, we get
$ (- \sin \theta ) \dfrac{d \theta} {dt} = (\omega_2 - \omega_1 ) \cdot (u_1 \times u_2) $
And this expression shows explicitly that even when the initial angle $\theta = 90^\circ$, then $\dfrac{d \theta}{dt} $ will be different from $0$ (as long as $\omega_2 - \omega_1 \ne 0 $).
If, on the other hand $u_1$ and $u_2$ were parallel initially, then the above equation cannot be used to calculate $\dfrac{d \theta}{dt} $, because we have $0$ on both sides of the equation ( $\sin \theta = 0 $ and $ u_1 \times u_2 = 0$ ).
In this case, we know that $u_1(t)$ and $u_2(t)$ will move on circles defined parametrcially by
$ u_1(t) = u \cos(\omega_1 t) + (\hat{\omega_1} \times u ) \sin(\omega_1 t) $
$ u_2(t) = u \cos(\omega_2 t) + ( \hat{\omega_2} \times u ) \sin(\omega_2 t ) $
Therefore
$ \cos(\theta(t) ) = u_1 \cdot u_2 = \cos(\omega_1 t) \cos(\omega_2 t) + (\hat{\omega_1} \times u ) \cdot ( \hat{\omega_2} \times u ) \sin(\omega_1 t ) \sin(\omega_2 t) $
Now,
$ (a \times c) \cdot (b \times c ) = c \cdot ( (a \times c) \times b )\\ = - c \cdot ( b \times (a \times c) ) = - c \cdot ( a (b.c) - c (a.b) )\\ = c^2 (a.b) - (a.c)(b.c) $
Therefore, $ (\hat{\omega_1} \times u ) \cdot ( \hat{\omega_2} \times u ) = \hat{\omega_1} \cdot \hat{\omega_2} $
And we now have
$ \cos(\theta(t)) = \cos(\omega_1 t) \cos(\omega_2 t) + (\hat{\omega_1} \cdot \hat{\omega_2}) \sin(\omega_1 t) \sin(\omega_2 t) $
The above equations models the angle $\theta(t)$ as a function of time $t$.
Using Taylor series expansion, you can show that
$ | \dfrac{ d\theta}{dt } | = \| \omega_1 - \omega_2 \| = \sqrt{ (\omega_1 - \omega_2 ) \cdot (\omega_1 - \omega_2 ) } $
The above expression gives the magnitude (the absolute value) of $\dfrac{d \theta}{dt} $ at $t = 0$ in the special case that initially the two vectors $u_1 $ and $u_2$ coincide.