If you forget about the angular velocity/angular momentum rigid body interpretation, the actual trajectories of the motion of the points are just circles on the sphere, traversed by the endpoints of the rod, which is a sphere of radius $l$ and centered at the center of mass $G$. The circles are simply the ones obtained by intersecting the sphere with planes. The motion of the point along such circle is uniform. This fact follows from the fact that the motion on the sphere is a curve under an acceleration always tangent to the sphere, always of the same magnitude and always perpendicular to the direction of motion. And the only curves that satisfy this condition are the circles described above, traversed uniformly.
But, if you still want to go over the angular velocity interpretation, with details, here it is.
First, I think that this situation is not equivalent to a spinning rod in a constant gravitational field.
Since the situation described in this problem is an abstract one, i.e. the rod is massless, connecting point masses and not rigid bodies, to treat it properly as a rigid body problem, first we need write the model for the rigid body version which is the problem for a pair of solid spheres of radius $r$, representing the mass-points, connected by a massless rod. Then, if you want, you can set $r = 0$ and get the solution for mass points.
It is a good idea to look at this problem from the coordinate system $G \, \vec{e}_1, \, \vec{e}_2, \, \vec{e}_3$ firmly attached to the rigid body (and thus moving with respect to an outside, world inertial coordinate system). The coordinate system $G \, \vec{e}_1, \, \vec{e}_2, \, \vec{e}_3$ is non-inertial. Let us assume that the origin of this coordinate system is at the center of mass $G$, the unit vector $\vec{e}_1$ is aligned with the segment from $G$ to one of the mass-points, while the other two unit vectors $\vec{e}_2$ and $\vec{e}_3$ are perpendicular to $\vec{e}_1$ and perpendicular to each other.
Spheres have the same moment of inertia $\frac{2}{5}mr^2$ with respect to any axis. Since inertia moments are additive, the inertia moment with respect to the $G\,\vec{e}_1$ axis is $I_1 = \frac{4}{5}mr^2$.
By the parallel axis theorem, the inertia moments with respect to the other two axes $G\,\vec{e}_2$ and $G\,\vec{e}_3$ are the same and equal to $I = I_1 + 2ml^2$
In other words, the inertia matrix in this body-fixed coordinate system takes the form
$$J = I_1 \, P_1 \, + I \, P_{23}$$ where $P_1$ is the orthogonal projection operator from 3D space onto the $G\vec{e}_1$ axis, and $P_{23}$ is the orthogonal projection matrix that maps the 3D space onto the coordinate plane spanned by $\vec{e}_2$ and $\vec{e}_3$, i.e.
$$P_{1} \vec{e}_1 = \vec{e}_1 \,\,\,\,\, P_{1} \vec{e}_2 = 0, \,\,\,\,\, P_{23} \vec{e}_3 = 0$$
$$P_{23} \vec{e}_1 = 0, \,\,\,\,\, P_{23} \vec{e}_2 = \vec{e}_2, \,\,\,\,\, P_{23} \vec{e}_3 = \vec{e}_3$$
In this coordinate system, the velocity vector is given by
$$\vec{v} = \vec{\omega} \times (\,l\,\vec{e}_1\,) = l\, \vec{\omega} \times \vec{e}_1$$
Hence, a unit vector, perpendicular to both $l \, \vec{e}_1$ and $\vec{v}$ is
$$\frac{ \big(\,\vec{\omega} \times \vec{e}_1 \,\big) \times \vec{e}_1}{|\vec{\omega} \times \vec{e}_1|}$$
And therefore, the two forces of constant magnitude $F$, as drawn on the picture, are
$$\vec{F} \, = \, F \, \frac{ \big(\,\vec{\omega} \times \vec{e}_1 \,\big) \times \vec{e}_1}{|\vec{\omega} \times \vec{e}_1|} \,\,\, \text{ and } \,\,\, -\, \vec{F} \, = \, - \, F \, \frac{ \big(\,\vec{\omega} \times \vec{e}_1 \,\big)\times \vec{e}_1}{|\vec{\omega} \times \vec{e}_1|}$$
The forces are applied at the points with position vectors $\vec{r} = l \, \vec{e}_1$ and $- \, \vec{r} = -\, l \, \vec{e}_1$ respectively. Consequently, the torques should be
\begin{align}
\vec{r} \times \vec{F} \, =&\, l \, \vec{e}_1 \times \left(\, F \, \frac{ \big(\, \vec{\omega} \times \vec{e}_1\,\big) \times \vec{e}_1}{|\vec{\omega} \times \vec{e}_1|} \,\right) \, = \, l F \,\, \frac{\vec{e}_1 \times \Big( \vec{e}_1 \times \big(\, \vec{e}_1 \times \vec{\omega}\,\big)\Big)}{| \vec{e}_1 \times \vec{\omega}|} \\
=& \,-\, lF \,\, \frac{ \,\vec{e}_1 \times \vec{\omega}\,}{|\vec{e}_1 \times \vec{\omega}|}
\end{align}
and
$$(-\vec{r}) \times (- \vec{F}) \, =\, \vec{r} \times \vec{F} \, = \, \,-\, lF \,\, \frac{\, \vec{e}_1 \times \vec{\omega}\,}{|\vec{e}_1 \times \vec{\omega}|}$$ respectively, so the total torque acting on the body is
$$-\, 2lF \,\, \frac{\, \vec{e}_1 \times \vec{\omega}\,}{|\vec{e}_1 \times \vec{\omega}|}$$
The equations of motion in the reference frame $G \, \vec{e}_1, \, \vec{e}_2, \, \vec{e}_3$ attached to the body, as described, are
\begin{align}
&\frac{d}{dt} J \, \vec{\omega} \, = \, \big(\,J \, \vec{\omega}\,\big) \times \vec{\omega} \, -\, 2lF \,\, \frac{\, \vec{e}_1 \times \vec{\omega}\,}{|\vec{e}_1 \times \vec{\omega}|}\\
&\frac{d}{dt} U^T \, = \, - \big(\,\vec{\omega} \times \cdot \,\big)\, U^T
\end{align}
where $\vec{\omega}$ is the angular velocity vector in the body-fixed coordinate frame $G \, \vec{e}_1, \, \vec{e}_2, \, \vec{e}_3$ and the matrix $U \, \in \, \text{SO}(3)$ is the rotation matrix that takes any coordinates from the body-fixed frame $G \, \vec{e}_1, \, \vec{e}_2, \, \vec{e}_3$ to the corresponding coordinates in the inertial world coordinate system, the latter also centered at the center of mass $G$. The notation $\big(\,\vec{\omega} \times \cdot \,\big)$ is the matrix that satisfies $$\big(\,\vec{\omega} \times \cdot \,\big) \,\vec{r} \, = \, \vec{\omega} \times \vec{r}$$ for any three dimensional vector $\vec{r}$ in the body-fixed frame.
Observe that the first (vector) equation is decoupled from the second (matrix) equation, so we will solve it separately.
First, let us set $\, \omega_1 \, \vec{e}_1 = P_1\, \vec{\omega}\,$ and $\,\vec{\omega}_{23} = P_{23}\, \vec{\omega}\,$. Then
$$\vec{\omega} = {\omega}_1 \, \vec{e}_1 \, + \, \vec{\omega}_{23}$$
Then
$$J \, \vec{\omega} \,=\,I_1 \, P_1 \vec{\omega} \, + \, I P_{23} \, \vec{\omega} = I_1 \, \omega_1\, \vec{e}_1 \, + \, I\, \vec{\omega}_{23}$$
and the equation becomes
\begin{align}
I_1 \, \frac{d\omega_1}{dt} \vec{e}_1 \, + \, I\, \frac{d}{dt}\vec{\omega}_{23} \, =& \,\big(\, I_1 \, \omega_1\, \vec{e}_1 \, + \, I\, \vec{\omega}_{23}\,\big) \times \big(\, \omega_1 \, \vec{e}_1 + \vec{\omega}_{23} \,\big)\\ & \, -\, 2lF \,\, \frac{\, \vec{e}_1 \times\big( {\omega}_1 \vec{e}_1 + \vec{\omega}_{23}\big)\,}{|\vec{e}_1 \times \big( {\omega}_1 \vec{e}_1 + \vec{\omega}_{23}\big)|}\\
=& \, \big(\, I_1 \, \omega_1\, \, - \, I\, \omega_1\,\big) \,\vec{e}_1 \times \vec{\omega}_{23} \, -\, 2lF \,\, \frac{\, \vec{e}_1 \times \vec{\omega}_{23}\,}{|\vec{e}_1 \times \vec{\omega}_{23}|}
\end{align}
or equivalently,
\begin{align}
I_1 \, \frac{d\omega_1}{dt} \vec{e}_1 \, + \, I\, \frac{d}{dt}\vec{\omega}_{23} \,
=& \, \left(\big(\, I_1 - I\,\big) \, \omega_1 \, - \,\, \frac{\, 2lF \,}{|\vec{e}_1 \times \vec{\omega}_{23}|}\right) \vec{e}_1 \times \vec{\omega}_{23}
\end{align}
From before, $I = I_1 + 2ml^2$, so $I - I_1 = 2ml^2$, which makes the equation
\begin{align}
I_1 \, \frac{d\omega_1}{dt} \vec{e}_1 \, + \, I\, \frac{d}{dt}\vec{\omega}_{23} \,
=& \, - \, \left( 2ml^2\, \omega_1 \, + \,\, \frac{\, 2lF \,}{|\vec{e}_1 \times \vec{\omega}_{23}|}\right) \vec{e}_1 \times \vec{\omega}_{23}
\end{align}
Next, let us simplify the norm $|\vec{e}_1 \times \vec{\omega}_{23}|$. By definition, $\vec{\omega}_{23}$ is the orthogonal projection of the angular velocity $\vec{\omega}$ onto the coordinate plane spanned by the coordinate vectors $\vec{e}_2$ and $\vec{e}_3$, which are perpendicular to the remaining coordinate vector $\vec{e}_1$. Therefore, $\vec{e}_1$ is perpendicular to $\vec{\omega}_{23}$ which is equivalent to the dot product $\vec{e}_1 \cdot \vec{\omega}_{23} \, = \, 0$. Hence the norm is
\begin{align}
|\vec{e}_1 \times \vec{\omega}_{23}| &= |\vec{e}_1| \, |\vec{\omega}_{23}| = |\vec{\omega}_{23}|
\end{align}
and the equation can be written in the from
\begin{align}
I_1 \, \frac{d\omega_1}{dt} \vec{e}_1 \, + \, I\, \frac{d}{dt}\vec{\omega}_{23} \,=& \, - \, \left( 2ml^2\, \omega_1 \, + \,\, \frac{\, 2lF \,}{|\vec{\omega}_{23}|}\right) \vec{e}_1 \times \vec{\omega}_{23}
\end{align}
Now, since the equation has two orthogonal components, one along $\vec{e}_1$ and the other parallel to $G \,\vec{e}_2 \, \vec{e}_3$, let us dot multiply both sides of the equation with $\vec{e}_1$
\begin{align}
\vec{e}_1 \cdot \left( I_1 \, \frac{d\omega_1}{dt} \vec{e}_1 \, + \, I\, \frac{d}{dt}\vec{\omega}_{23}\right) \,=& \, - \, \vec{e}_1 \cdot \left( \, \left( 2ml^2\, \omega_1 \, + \,\, \frac{\, 2lF \,}{|\vec{\omega}_{23}|}\right) \vec{e}_1 \times \vec{\omega}_{23} \right)\\
&\\
I_1 \, \frac{d\omega_1}{dt} (\vec{e}_1 \cdot \vec{e}_1) \, + \, I\, \left(\vec{e}_1 \cdot \frac{d}{dt}\vec{\omega}_{23} \right) \,=& \, - \, \left( 2ml^2\, \omega_1 \, + \,\, \frac{\, 2lF \,}{|\vec{\omega}_{23}|}\right) \big(\, \vec{e}_1 \cdot ( \,\vec{e}_1 \times \vec{\omega}_{23} ) \big)\\
\end{align}
However, by the properties of the cross-product $\vec{e}_1$ is perpendicular to $\vec{e}_1 \times \vec{\omega}_{23}$ which means that $$\big(\, \vec{e}_1 \cdot ( \,\vec{e}_1 \times \vec{\omega}_{23} ) \big) = 0$$ Consequently, we arrive at
$$I_1 \, \frac{d\omega_1}{dt} (\vec{e}_1 \cdot \vec{e}_1) \, + \, I\, \frac{d}{dt}(\vec{e}_1 \cdot \vec{\omega}_{23} ) \,= \, 0$$
and since $(\vec{e}_1 \cdot \vec{e}_1) = 1$ and $(\vec{e}_1 \cdot \vec{\omega}_{23} ) = 0$ we conclude that
$$I_1 \, \frac{d\omega_1}{dt} = 0$$
Since we have assumed that $I_1$ is not zero, the only remaining possibility is that
$$\frac{d\omega_1}{dt} = 0$$
which means that the $\vec{e}_1$ component $\omega_1$ of the angular velocity is a constant. Let us denote this constant by $\omega_1^0$
The vector differential equation for the angular velocity reduces to
\begin{align}
I\, \frac{d}{dt}\vec{\omega}_{23} \,=& \, - \, \left( 2ml^2\, \omega_1^0 \, + \,\, \frac{\, 2lF \,}{|\vec{\omega}_{23}|}\right) \vec{e}_1 \times \vec{\omega}_{23}
\end{align}
The final simplification of the equation comes from deducing that the magnitude of $\vec{\omega}_{23}$ is constant. To prove that, we begin by differentiating with respect to time $t$ the magnitude of $\vec{\omega}_{23}$ squared:
\begin{align}
\frac{d}{dt} \left(\frac{1}{2} \,|\vec{\omega}_{23}|^2\right) \, =& \, \vec{\omega}_{23} \cdot \frac{d\vec{\omega}_{23}}{dt}\\
=& \, \vec{\omega}_{23} \cdot \left( - \, \frac{1}{I}\left( 2ml^2\, \omega_1^0 \, + \,\, \frac{\, 2lF \,}{|\vec{\omega}_{23}|}\right) \vec{e}_1 \times \vec{\omega}_{23} \right)\\
=& \, - \, \frac{1}{I}\left( 2ml^2\, \omega_1^0 \, + \,\, \frac{\, 2lF \,}{|\vec{\omega}_{23}|}\right) \big( \vec{\omega}_{23} \cdot ( \vec{e}_1 \times \vec{\omega}_{23}) \big)\\
=& \, 0
\end{align}
because by the properties of the cross-product $\vec{\omega}_{23}$ is perpendicular to $\vec{e}_1 \times \vec{\omega}_{23}$. the latter chain of identities is possible only when
$$\frac{1}{2} \,|\vec{\omega}_{23}|^2 = \text{const}$$
which means that the magnitude $|\vec{\omega}_{23}|$ is also a constant. Let us denote this constant by
$$\omega_{23}^0 = |\vec{\omega}_{23}|$$
Finally, the vector differential equation for the angular velocity simplifies to the linear constant coefficient vector differential equation
\begin{align}
\frac{d}{dt}\vec{\omega}_{23} \,=& \, - \, \left( \frac{2ml^2}{I}\, \omega_1^0 \, + \,\, \frac{\, 2lF \,}{I \,{\omega}_{23}^0}\right) \vec{e}_1 \times \vec{\omega}_{23}
\end{align}
All together, the total system reduces to
\begin{align}
&\vec{\omega} = {\omega_1^0} \,\vec{e}_1 + \,\vec{\omega}_{23}\\
&\omega_{23}^0 = |\vec{\omega}_{23}|\\
&I = \frac{4}{5}mr^2 \, + \, 2ml^2\\
&\nu_0 = - \, \left( \frac{2ml^2}{I}\, \omega_1^0 \, + \,\, \frac{\, 2lF \,}{I \,{\omega}_{23}^0}\right)\\
&\frac{d}{dt}\vec{\omega}_{23} \,= \, \nu_0\, \vec{e}_1 \times \vec{\omega}_{23}\\
&\frac{d}{dt} U^T \, = \, - \big(\,\vec{\omega} \times \cdot \,\big)\, U^T
\end{align}
At this point, we can solve the vector differential equation for $\vec{\omega}_{12}$, where the solution is
$$\vec{\omega}_{12} \, = \, \omega_{23}^0\cos(\nu_0 t)\,\vec{e}_2 \, + \, \omega_{23}^0\sin(\nu_0 t)\,\vec{e}_3 $$
and so the full angular velocity is
$$\vec{\omega} \, = \, \omega_1^0\, \vec{e}_1 \, + \, \omega_{23}^0\cos(\nu_0 t)\,\vec{e}_2 \, + \, \omega_{23}^0\sin(\nu_0 t)\,\vec{e}_3 $$
In your case, $\omega_1^0 = 0$ (and if you wish, you can set $r=0$, i.e. shrink the pair of spheres to mass points) the angular velocity is
$$\vec{\omega} \, = \, \omega_0\cos(\nu_0 t)\,\vec{e}_2 \, + \, \omega_0\sin(\nu_0 t)\,\vec{e}_3 $$
where $\omega_0$ is the constant magnitude of the angular velocity. As you can see, $\vec{\omega}$ traverses a circle, located in the coordinate plane $\,G \, \vec{e}_2 \, \vec{e}_3\,$, centered at $G$, and of radius $\omega_0$. Moreover, $\vec{\omega}$ traverses that circle uniformly, at a constant frequency of $\nu_0$.
Let us consider a rigid body consisting of two particles of mass $m$ at
either end of a rigid massless rod of length $2l$.
