I don't know if this question still matters, but I will answer it. Basically, the idea is very general and one can derive all the equations of motion from first principles.
Given a system of mass points with mass $m_j \, :\, j =1 ... n$ and position vectors $\vec{r}_j \, :\, j =1 ... n$ in an inertial coordinate system they satisfy Newton's equations of mottion
$$m_j \, \frac{d^2\vec{r}_j}{dt^2} = \vec{f}_j\Big(\, \vec{r}_1, ..., \vec{r}_n, \frac{d\vec{r}_1}{dt}, ..., \frac{d\vec{r}_n}{dt}, t\, \Big) \text{ for } j = 1...n$$
where $$\vec{f}_j = \vec{f}_j\Big(\, \vec{r}_1, ..., \vec{r}_n, \frac{d\vec{r}_1}{dt}, ..., \frac{d\vec{r}_n}{dt}, t\, \Big)$$ are the forces acting on each particle. We can cross-product multiply both sides of each equation as
$$m_j\left( \vec{r}_j \times \frac{d^2\vec{r}_j}{dt^2} \right) = \vec{r}_j \times \vec{f}_j$$ Due to the properties of the cross-product
$$\frac{d}{dt}\left( \vec{r}_j \times \frac{d\vec{r}_j}{dt} \right) \, =\, \frac{d\vec{r}_j}{dt} \times \frac{d\vec{r}_j}{dt} \, + \, \vec{r}_j \times \frac{d^2\vec{r}_j}{dt^2}\, = \, \vec{r}_j \times \frac{d^2\vec{r}_j}{dt^2} $$ Thus, we can rewrite the equations above as
$$ \frac{d}{dt}\, \left(\, m_j \Big( \vec{r}_j \times \frac{d\vec{r}_j}{dt} \Big)\, \right) \, =\, \vec{r}_j \times \vec{f}_j $$ Finally, we can sum them together to obtain
$$ \sum_{j=1}^{n} \, \frac{d}{dt}\, \left(\, m_j \Big( \vec{r}_j \times \frac{d\vec{r}_j}{dt} \Big)\, \right) \, =\, \sum_{j=1}^{n} \, \vec{r}_j \times \vec{f}_j $$
Now, let us focus on the rotating bar. Fix an inertial coordinate system $O\,\vec{e}_x \,\vec{e}_y\,\vec{e}_z $ with origin $O$ coinciding with the point of rotation of the bar. We represent the rotating rod as a system of continuum many particles, each of which is represented by the position vector $\vec{r} = \vec{r}(t) = x(t)\, \vec{e}_x + y(t)\, \vec{e}_y + z(t)\, \vec{e}_z $ pointing from $O$ to the point on the bar which represents the particles at time $t$. By $\mu(\vec{r})$ we represent the mass distribution (mass-density) of the bar. The important part is that the bar is a rigid body, so we can consider a coordinate system $O\,\vec{E}_X \,\vec{E}_Y\,\vec{E}_Z$ attached firmly to the bar, which means it rotates together with the bar and the bar is at rest with respect to $O\,\vec{E}_X \,\vec{E}_Y\,\vec{E}_Z$. The position of a point on the bar can be expressed as $\vec{R} = X\, \vec{E}_X + Y\, \vec{E}_Y + Z\, \vec{E}_Z $ wher $\vec{R}$ doesn't change with time with respect to $O\,\vec{E}_X \,\vec{E}_Y\,\vec{E}_Z$. Then, there is a time dependent rotation matrix $U = U(t)\, \in \, \text{SO}(3)$ such that
$$\vec{r}(t) = U(t)\,\vec{R}$$ where
$$\vec{r}(t) = \begin{bmatrix} x(t)\\y(t)\\z(t)\end{bmatrix} \, \text{ and } \, \vec{R} = \begin{bmatrix} X\\Y\\Z\end{bmatrix}$$
Thus, as already explained before, for each point from the rigid bar in the inertial coordinate system $O\,\vec{e}_x \,\vec{e}_y\,\vec{e}_z$ we have $$\frac{d}{dt}\, \left(\, \mu(\vec{r}) \Big( \vec{r} \times \frac{d\vec{r}}{dt} \Big)\, \right) \, =\, \vec{r} \times \vec{f}\Big(\vec{r}, \frac{d\vec{r}}{dt}, t\Big) $$ where
$\vec{f} = \vec{f}\Big(\vec{r}, \frac{d\vec{r}}{dt}, t\Big)$ is the force acting on the bar at the point $\vec{r}$. The idea is to sum all the points on the bar up, but it is difficult because they change position with time. Therefore, switching to the system $O\,\vec{E}_X \,\vec{E}_Y\,\vec{E}_Z$ the position vectors $\vec{R}$ do not change with time. Thus, we calculate
$$\frac{d\vec{r}}{dt} = \frac{d U}{dt}\, \vec{R} = U \big(\vec{\Omega} \times \vec{R}\big)$$ because for any time-dependent orthogonal matrix $U = U(t)$, there exists a time-dependent vector $\vec{\Omega} = \vec{\Omega}(t)$, called angular velocity, such that $$U^{-1}\frac{d U}{dt} \vec{R} = U^T\frac{d U}{dt} \vec{R} = \vec{\Omega} \times \vec{R}$$ Consequently,
\begin{align}\frac{d}{dt}\, \left(\, \mu(\vec{r}) \Big( \vec{r} \times \frac{d\vec{r}}{dt} \Big)\, \right) \, =& \, \frac{d}{dt}\, \left(\, \mu(\vec{R}) \Big( U\,\vec{R} \times U\, \big(\vec{\Omega} \times \vec{R}\big)\,\Big)\, \right) = \mu(\vec{R})\, \frac{d}{dt}\, \left(\, U \,\Big(\vec{R} \times \big(\vec{\Omega} \times \vec{R}\big)\,\Big)\, \right) \end{align}
To simplify the notation, let us set the following linear transformation acting linearly on $\vec{\Omega}$ and changing quadratically with respect to $\vec{R}$
$$A\big(\vec{R}\big)\, \vec{\Omega} = \vec{R} \times \big(\vec{\Omega} \times \vec{R}\big) = |\vec{R}|^2\, \vec{\Omega} - \big(\vec{R} \circ \vec{\Omega}\big)\, \vec{R}$$
However, recall that $\vec{R}$ does not change with time $t$. Thus
\begin{align}
\frac{d}{dt}\, \left(\, \mu(\vec{r}) \Big( \vec{r} \times \frac{d\vec{r}}{dt} \Big)\, \right) \, =& \, \mu(\vec{R})\, \frac{d}{dt}\, \left(\, U \,\Big(\vec{R} \times \big(\vec{\Omega} \times \vec{R}\big)\,\Big)\, \right) \, = \, \mu(\vec{R})\, \frac{d}{dt}\, \left(\, U \, A\big(\vec{R}\big)\, \vec{\Omega}\,\right)\\
=& \,\mu(\vec{R})\, \left(\, \frac{dU}{dt} \, A\big(\vec{R}\big)\, \vec{\Omega} \, + \, U \, A\big(\vec{R}\big)\, \frac{d\vec{\Omega}}{dt} \, \right)\\
=& \,\mu(\vec{R})\, \left(\, U \,\Big[ \vec{\Omega} \times \Big(\, A\big(\vec{R}\big)\, \,\vec{\Omega}\,\Big)\,\Big]\, + \, U \, A\big(\vec{R}\big)\, \frac{d\vec{\Omega}}{dt} \, \right)\\
=& \, U\, \left(\, \Big[ \vec{\Omega} \times \Big(\,\mu(\vec{R})\, A\big(\vec{R}\big)\, \,\vec{\Omega}\,\Big)\,\Big]\, + \, \mu(\vec{R})\, A\big(\vec{R}\big)\, \frac{d\vec{\Omega}}{dt} \, \right)\\
=&\, U\, \left(\,\mu(\vec{R})\, A\big(\vec{R}\big)\, \frac{d\vec{\Omega}}{dt}\, + \,\Big[ \vec{\Omega} \times \Big(\,\mu(\vec{R})\, A\big(\vec{R}\big)\, \,\vec{\Omega}\,\Big)\,\Big] \,\right)
\end{align} On the other hand, the torque is
\begin{align}\vec{r} \times \vec{f}\Big(\vec{r}, \frac{d\vec{r}}{dt}, t \Big) \, =& \,
\Big[U\,\vec{R}\Big] \times \Big[ \, U\, U^{-1}\,\vec{f}\Big(\, U\vec{R}, \, U\, \big(\vec{\Omega}\times \vec{R}\big), \, t \,\Big)\, \Big] \\ =& U \Big[\, \vec{R} \times U^{-1}\,\vec{f}\Big(\, U\vec{R}, \, U\, \big(\vec{\Omega}\times \vec{R}\big), \, t \,\Big) \, \Big]\\
=& \, U \Big[\, \vec{R} \times U^T\vec{f}\Big(\, U\vec{R}, \, U\, \big(\vec{\Omega}\times \vec{R}\big), \, t \,\Big) \, \Big]\\
=& \, U \Big[\, \vec{R} \times \vec{F}\Big(\, \vec{R},\, U,\, \, \vec{\Omega}, \, t \,\Big) \, \Big]
\end{align} Where $ \vec{F}\Big(\, \vec{R},\, U,\, \, \vec{\Omega}, \, t \,\Big) = U^T\vec{f}\Big(\, U\vec{R}, \, U\, \big(\vec{\Omega}\times \vec{R}\big), \, t \,\Big)$. Consequently, the equations of motion
$$\frac{d}{dt}\, \left(\, \mu(\vec{r}) \Big( \vec{r} \times \frac{d\vec{r}}{dt} \Big)\, \right) \, = \, \vec{r} \times \vec{f}\Big(\vec{r}, \frac{d\vec{r}}{dt}, t \Big)$$ can be written as
$$ U\, \left(\,\mu(\vec{R})\, A\big(\vec{R}\big)\, \frac{d\vec{\Omega}}{dt}\, + \,\Big[ \vec{\Omega} \times \Big(\,\mu(\vec{R})\, A\big(\vec{R}\big)\, \,\vec{\Omega}\,\Big)\,\Big] \,\right) \, = \, U \Big[\, \vec{R} \times \vec{F}\Big(\, \vec{R},\, U,\, \, \vec{\Omega}, \, t \,\Big) \, \Big]$$ When we multiply both sides of the equations with the inverse orthogonal matrix $U^{-1} = U^T$ we obtain the equations written in the rotating frame $O\,\vec{E}_X \vec{E}_Y \vec{E}_Z$, firmly attached to the rigid bar
$$ \mu(\vec{R})\, A\big(\vec{R}\big)\, \frac{d\vec{\Omega}}{dt}\, + \,\Big[ \vec{\Omega} \times \Big(\,\mu(\vec{R})\, A\big(\vec{R}\big)\, \,\vec{\Omega}\,\Big)\,\Big] \, = \, \vec{R} \times \vec{F}\Big(\, \vec{R},\, U,\, \, \vec{\Omega}, \, t \,\Big) $$
However, these are the equations of motion for only one point $\vec{R}$ from the bar. Summing them up means to integrate them, obtaining
$$ \int_{B}\, \left(\, \mu(\vec{R})\, A\big(\vec{R}\big)\, \frac{d\vec{\Omega}}{dt}\, + \,\Big[ \vec{\Omega} \times \Big(\,\mu(\vec{R})\, A\big(\vec{R}\big)\, \,\vec{\Omega}\,\Big)\,\Big] \, \right) \, dR \, = \, \int_{B}\, \left(\, \vec{R} \times \vec{F}\Big(\, \vec{R},\, U,\, \, \vec{\Omega}, \, t \,\Big)\, \right) \, dR $$ where $dR$ is the Lebesgue measure on the rigid body we are dealing with. The linearity of the integral yields
$$ \int_{B}\, \left(\, \mu(\vec{R})\, A\big(\vec{R}\big)\, \frac{d\vec{\Omega}}{dt}\, \right) \, dR \, + \, \vec{\Omega} \times \int_{B}\, \left(\,\mu(\vec{R})\, A\big(\vec{R}\big)\, \,\vec{\Omega}\,\Big) \, \right) \, dR \, = \, \int_{B}\, \left(\, \vec{R} \times \vec{F}\Big(\, \vec{R},\, U,\, \, \vec{\Omega}, \, t \,\Big)\, \right) \, dR $$ Then we can define the linear transformation (linear matrix) $J$ as follows
$$J\, \vec{\Omega} = \int_{B}\, \left(\,\mu(\vec{R})\, A\big(\vec{R}\big)\, \,\vec{\Omega}\,\Big) \, \right) \, dR$$
$$J\,\frac{d\vec{\Omega}}{dt} = \int_{B}\, \left(\, \mu(\vec{R})\, A\big(\vec{R}\big)\, \frac{d\vec{\Omega}}{dt}\, \right) \, dR $$ and the expression for the total torque is
$$\vec{T}\big(U, \vec{\Omega}, t\big) \, = \, \int_{B}\, \left(\, \vec{R} \times \vec{F}\Big(\, \vec{R},\, U,\, \, \vec{\Omega}, \, t \,\Big)\, \right) \, dR$$ Thus, we have arrived at the standard equations of motion for the rigid body, which are
\begin{align}
J \, &\frac{d\vec{\Omega}}{dt} + \vec{\Omega} \times J\, \vec{\Omega} \, = \,\vec{T}\big(U, \vec{\Omega}, t\big)\\
& \frac{dU}{dt} = U \, (\vec{\Omega} \times \cdot)
\end{align} where by $(\vec{\Omega} \times \cdot)$ I have denoted the skew-symmetric matrix that acts on any vector $\vec{R}$ as
$$(\vec{\Omega} \times \cdot) \,\vec{R} = \vec{\Omega} \times \vec{R}$$ The matrix $J$ is what is called the inertia tensor. These equations are very general equations of motion for any rigid body rotating around a fixed point $O$.
In your case, the body $B$ is a bar, so one dimensional segment of length $l$. We can assume that the coordinate system $O\,\vec{E}_X \vec{E}_Y \vec{E}_Z$ is attached to the bar so that the bar is positioned along the $X$ axis only, so any vector $\vec{R} = X \, \vec{E}_X$. Furthermore, the rotation happens so that the $z-$axis stays fixed which means that $\vec{e}_z = \vec{E}_Z$ is fixed. The rotations that keep the $z-$axis fixed have the form
$$U = U(\theta) = \begin{bmatrix} \cos(\theta) & - \, \sin(\theta) & 0 \\
\sin(\theta) & \cos(\theta) & 0 \\
0 & 0 & 1
\end{bmatrix}$$ where $\theta = \theta(t)$ determines the matrix's change with respect to time, and consequently the angular velocity is
$$\vec{\Omega}\times \vec{R} = \Big(U(\theta)^T \frac{d}{dt}U(\theta) \Big)\, \vec{R} = \frac{d\theta}{dt}\, \vec{E}_Z \times \vec{R}$$ for any vector $\vec{R}$, i.e. $\vec{\Omega} = \frac{d\theta}{dt}\, \vec{E}_Z = \omega\, \vec{E}_Z $ where $\omega = \frac{d\theta}{dt}$.
First, let us calculate the torques. If you have to account for gravity, say pointing down the $y-$axis of the inertial coordinate system $O\, \vec{e}_x\vec{e}_y\vec{e}_z$, then the gravitational force is
$$\vec{f}_{gr} = -\, mg\,\vec{e}_y$$ and transforms to the rotating system $O\, \vec{E}_X\vec{E}_Y\vec{E}_Z$ as
$$\vec{F}_{gr} = U^{-1}\, \vec{f}_{gr} = U^T\, \vec{f}_{gr} = -\, mg\,U^T\,\vec{e}_y = -\, mg\, \big(\sin(\theta)\, \vec{E}_X + \cos(\theta)\, \vec{E}_Y\big)$$ because
$$\vec{e}_y =\begin{bmatrix} 0\\ 1 \\0 \end{bmatrix} \text{ and } \, U^T\,\vec{e}_y = \begin{bmatrix} \cos(\theta) & - \, \sin(\theta) & 0 \\
\sin(\theta) & \cos(\theta) & 0 \\
0 & 0 & 1
\end{bmatrix}^T \begin{bmatrix} 0\\ 1 \\0 \end{bmatrix} = \begin{bmatrix} \cos(\theta) & \sin(\theta) & 0 \\
-\, \sin(\theta) & \cos(\theta) & 0 \\
0 & 0 & 1
\end{bmatrix} \begin{bmatrix} 0\\ 1 \\0 \end{bmatrix} = \begin{bmatrix} \sin(\theta)\\ \cos(\theta) \\0 \end{bmatrix}$$ Thus
$$\vec{R} \times \vec{F}_{gr} = - mg\, X \, \vec{E}_X \times \big(\sin(\theta)\, \vec{E}_X + \cos(\theta)\, \vec{E}_Y\big) = - \, mg\, X \, \cos(\theta) \, \vec{E}_X \times \vec{E}_Y = - \, mg\, X \, \cos(\theta) \, \vec{E}_Z$$
Thus, the total gravitational torque is
$$\vec{T}_{gr} = - \, mg\, \int_{0}^{l} \, \Big(\, X \,\cos(\theta) \, \vec{E}_Z \,\Big)dX = - \, mg\, \Big( \int_{0}^{l} \, X \, dX \Big)\, \cos(\theta) \, \vec{E}_Z = -\, \frac{mgl^2}{2}\, \cos(\theta) \, \vec{E}_Z$$ Assume, the fluid resistance force is modeled in the inertial coordinate system $O\, \vec{e}_x\vec{e}_y\vec{e}_z$ as
$$\vec{f} = - k \, \left|\frac{d\vec{r}}{dt}\right|^{\beta}\,\frac{d\vec{r}}{dt} $$
Having in mind that
$$\frac{d\vec{r}}{dt} = U \big(\vec{\Omega} \times \vec{R}\big) = U \Big(\big[\omega\, \vec{E}_Z\big] \times \big[X\,\vec{E}_X\big]\Big) = \omega\,X \, U \big(\vec{E}_Z \times \vec{E}_X\big) = \omega\,X \, U \,\vec{E}_Y $$
$$\vec{f} = - k \, \left|\frac{d\vec{r}}{dt}\right|^{\beta}\,\frac{d\vec{r}}{dt} =
- k \, \left|\omega\,X \, U \,\vec{E}_Y\right|^{\beta}\,\omega\,X \, U \,\vec{E}_Y =
- k \,U \left( \left|\omega\,X \,\vec{E}_Y\right|^{\beta}\,\omega\,X \,\vec{E}_Y\right) $$ so
$$\vec{F} = U^{-1}\vec{f} = - k \, \left|\omega\,X \,\vec{E}_Y\right|^{\beta}\,\omega\,X \,\vec{E}_Y = -\,k\, \omega^{\beta+1}\, X^{\beta+1} \big|\vec{E}_Y\big|^{\beta}\,\vec{E}_Y = -\,k\, \omega^{\beta+1}\, X^{\beta+1}\,\vec{E}_Y$$ Thus
$$\vec{R} \times \vec{F} = \big[ X \, \vec{E}_X\big] \times \big[ -\,k\, \omega^{\beta+1}\, X^{\beta+1}\,\vec{E}_Y \big] = -\,k\, \omega^{\beta+1} \, X^{\beta+2} \, \big(\vec{E}_X \times \vec{E}_Y\big) = -\,k\, \omega^{\beta+1} \, X^{\beta+2} \, \vec{E}_Z$$ Finally the total resistance torque is
$$\vec{T} = -\,k\, \int_{0}^{l}\Big(\, \omega^{\beta+1} \, X^{\beta+2} \, \vec{E}_Z \,\Big)dX = -\,k\, \Big(\int_{0}^{l} \, X^{\beta+2}dX \Big)\, \omega^{\beta+1}\, \vec{E}_Z = -\,\frac{k\,l^{\beta + 3}}{\beta+3}\, \omega^{\beta+1}\, \vec{E}_Z $$
To calculate the inertia tensor, one goes back to the expression
\begin{align}
A\big(\vec{R}\big)\,\vec{\Omega} =& |\vec{R}|^2\, \vec{\Omega} - \big(\vec{R} \circ \vec{\Omega}\big)\, \vec{R} = |X\, \vec{E}_X|^2\, \omega\, \vec{E}_Z - \big(X\,\vec{E}_X \circ \omega\,\vec{E}_Z\big)\, X\,\vec{E}_X\\
=& X^2\,\omega\,|\vec{E}_X|^2\, \vec{E}_Z - X\,\omega\,\big(\vec{E}_X \circ \vec{E}_Z\big)\, X\,\vec{E}_X\\ =& X^2\,\omega\, \vec{E}_Z
\end{align} because the $X-$ and $Z-$axes are orthogonal and therefore $\big(\vec{E}_X \circ \vec{E}_Z\big) = 0$ as well as $|\vec{E}_X|^2 = 1$. Analogously
$$A\big(\vec{R}\big)\,\frac{\vec{\Omega}}{dt} = X^2\,\frac{d\omega}{dt}\, \vec{E}_Z$$ Moreover, if we assume that the bar has length $l$ and mass $m$ which is distributed homogeneously along the bar, $\mu(\vec{R}) = \frac{m}{l}$.
Thus
$$ J\, \vec{\Omega} = \int_{B}\, \mu(\vec{R})\, A(\vec{R})\,\vec{\Omega}\, dR = \int_{0}^{l}\, \Big(\frac{m}{l}\, X^2\, dX \Big)\, {\omega}\, \vec{E}_Z = \frac{m}{l}\, \frac{l^3}{3}\, {\omega}\, \vec{E}_Z = \frac{m\,l^2}{3} \, \omega\, \vec{E}_Z$$
$$ J\, \frac{d\vec{\Omega}}{dt} = \int_{B}\, \mu(\vec{R})\, A(\vec{R})\,\frac{d\vec{\Omega}}{dt}\, dR = \int_{0}^{l}\, \Big(\frac{m}{l}\, X^2\, dX \Big)\, \frac{d\omega}{dt}\, \vec{E}_Z = \frac{m}{l}\, \frac{l^3}{3}\, \frac{d\omega}{dt}\, \vec{E}_Z = \frac{m\,l^2}{3} \, \frac{d\omega}{dt}\, \vec{E}_Z$$ Since in this case
$$\vec{\Omega} \times J\,\vec{\Omega} = \omega\, \vec{E}_Z \times \Big(\frac{m\,l^2}{3} \, {\omega}\, \vec{E}_Z \Big) = \Big(\frac{m\,l^2}{3}\, \omega^2\Big)\, \vec{E}_Z \times \vec{E}_Z =\vec{0}$$
the general equations of motion of the system
\begin{align}
J \, &\frac{d\vec{\Omega}}{dt} + \vec{\Omega} \times J\, \vec{\Omega} \, = \,\vec{T}\big(U, \vec{\Omega}, t\big)\\
& \frac{dU}{dt} = U \, (\vec{\Omega} \times \cdot)
\end{align} where by $(\vec{\Omega} \times \cdot)$ reduce to
\begin{align}
\frac{m\,l^2}{3} \, &\frac{d\omega}{dt}\, \vec{E}_Z \, = \, - \,\left(\frac{k\,l^{\beta+3}}{\beta+3}\right)\, \omega^{\beta+1}\, \vec{E}_Z \,-\,\frac{m\, g\, l^2}{2}\,\cos(\theta)\, \vec{E}_Z \\
& \frac{d\theta}{dt} = \omega
\end{align} which, after equating the coefficients in front of the vector $\vec{E}_Z$, yield
\begin{align}
\left(\frac{m\,l^2}{3}\right) \, &\frac{d\omega}{dt}\, = \, - \,\left(\frac{k\,l^{\beta+3}}{\beta+3}\right)\, \omega^{\beta+1}\,-\,\left(\frac{m\, g\, l^2}{2}\right)\,\cos(\theta) \\
& \frac{d\theta}{dt} = \omega
\end{align} or as one equation
\begin{align}
&\left(\frac{m\,l^2}{3}\right) \, \frac{d^2\theta}{dt^2}\, = \, - \,\left(\frac{k\,l^{\beta+3}}{\beta+3}\right)\, \left(\frac{d\theta}{dt}\right)^{\beta+1}\,-\,\left(\frac{m\, g\, l^2}{2}\right)\,\cos(\theta)\end{align}
Now, if the bar is horizontal and the gravity doesn't apply, then the equations are
\begin{align}
\left(\frac{m\,l^2}{3}\right) \, &\frac{d\omega}{dt}\, = \, - \,\left(\frac{k\,l^{\beta+3}}{\beta+3}\right)\, \omega^{\beta+1}\\
& \frac{d\theta}{dt} = \omega
\end{align} and if you want the bar to rotate at (or very near to) a fixed angular velocity $\omega_0 \, \vec{E}_Z$, you can choose the controlling torque $\vec{T}_{c}$ that cancels the fluid resistance to be say
$$\vec{T}_{c} = \left(\,- \, K \, \big(\omega - \omega_0\big) \, + \, \left( \frac{k\,l^{\beta+3}}{\beta+3}\right)\, \omega^{\beta+1}\, \right) \, \vec{E}_Z$$ where $K>0$. This torque also adds Lyapunov asymptotic stability around $\omega_0$
