Heavy Symmetrical Top

Alexei Gilchrist

The equations of motion for a heavy symmetrical top are derived in terms of the Euler angles and using a Lagrangian approach. The solutions are examined numerically.

1 Approach

Imagine that we are describing the motion of a heavy symmetrical top that is spinning around its axis of symmetry. There are potentially several ways we can tackle this problem, the approach we will take here is to derive the equations of motion in terms of Euler angles using the Lagrangian approach. The motivation of this approach is that we can consider two co-ordinate axes—one is the stationary axes we are observing from and in which we want to express our results. The other is the principal axes of the body for which we know the moments of inertia and can calculate quantities like the kinetic energy simply. The Euler-Lagrange equations will apply for any choice of generalised coordinates so this approach will work fine with the Euler angles once we can calculate the kinetic energy of the top.

2 Calculating \(\vec{\omega}\)

In the diagram below, an object with principal axes \(\{\hat{{x}}_1, \hat{{x}}_2, \hat{{x}}_3\}\) is rotated with respect to a fixed set of axes \(\{\hat{{x}}, \hat{{y}}, \hat{{z}}\}\) by the set of Euler angles \(\phi\), \(\theta\), and \(\psi\). The Euler angles can be thought of as a sequence of rotations—first around \(\hat{{x}}_3\) by \(\phi\), then around \(\hat{{x}}_1\) by \(\theta\), and finally around \(\hat{{x}}_3\) again by \(\psi\).

The angular velocity vectors for the Euler angles.

In order to calculate the kinetic energy of the object, we need to find the total angular velocity \(\vec{\omega}\) in the \(\{\hat{{x}}_1, \hat{{x}}_2, \hat{{x}}_3\}\) co-ordinate system (where the moment of inertia will be simple). Since we can think of the Euler angles as a nested sequence of rotations, we can add the angular velocities of each rotation,

\[\begin{equation*}\vec{\omega} = \dot{\vec{\phi}}+\dot{\vec{\theta}}+\dot{\vec{\psi}},\end{equation*}\]
and the task is to find the components of the angular velocities in the principal axes.

Each of the angular velocity vectors will be orthogonal to the plane of rotation with the direction given by the right-hand-rule. Clearly \(\vec{\dot{\psi}}\) points along \(\hat{{x}}_3\) so that

\[\begin{equation*}\dot{\vec{\psi}} = \begin{pmatrix} 0\\0\\ \dot{\psi} \end{pmatrix}.\end{equation*}\]
Similarly, \(\dot{\vec{\theta}}\) points along the join in the two planes depicted in the figure, and it will have no component in the \(\hat{{x}}_3\) direction. It’s components are split between \(\hat{{x}}_1\) and \(\hat{{x}}_2\). Projecting onto these axes we get,
\[\begin{equation*}\dot{\vec{\theta}} = \begin{pmatrix} \dot{\theta}\cos\psi \\ -\dot{\theta}\sin\psi \\ 0 \end{pmatrix}.\end{equation*}\]
Finally, \(\dot{\vec{\phi}}\) points along \(\hat{{z}}\) and it will have components along all the principal axes. The projection of \(\hat{{z}}\) to \(\hat{{x}}_3\) is \(\cos\theta\) which means that the projection onto the dashed line on the \(\psi\) plane of rotation in the figure will be \(\sin\theta\). This projection on the plane has components in the \(\hat{{x}}_1\) and \(\hat{{x}}_2\) direction given by \(\psi\) so that,
\[\begin{equation}\label{eq:z-axis} \hat{{z}} = \begin{pmatrix} \sin\theta\sin\psi \\ \sin\theta\cos\psi \\ \cos\theta \end{pmatrix},\end{equation}\]
and \(\dot{\vec{\phi}} = \dot{\phi}\hat{{z}}\).

Collecting results we get

\[\begin{equation}\label{eq:omega} \vec{\omega} = \begin{pmatrix} \dot{\phi}\sin\theta\sin\psi+\dot{\theta}\cos\psi \\ \dot{\phi}\sin\theta\cos\psi -\dot{\theta}\sin\psi \\ \dot{\phi}\cos\theta+\dot{\psi} \end{pmatrix}.\end{equation}\]
Using this result we can easily calculate the angular momentum since the inertia tensor is diagonal in the principal axes:
\[\begin{equation}\label{eq:L} \vec{L} = \mathbf{I}\vec{\omega} = \begin{pmatrix} I_1(\dot{\phi}\sin\theta\sin\psi+\dot{\theta}\cos\psi) \\ I_2(\dot{\phi}\sin\theta\cos\psi -\dot{\theta}\sin\psi) \\ I_3(\dot{\phi}\cos\theta+\dot{\psi}) \end{pmatrix}.\end{equation}\]

3 Equations of Motion

The top we are describing is symmetrical, meaning that along the moment’s of inertia along the principal axes are \(I_1=I_2 \ne I_3\), where \(\hat{{x}}_3\) is the axis of symmetry.

A symmetrical top spinning about its axis of symmetry.

Referring to the diagram above the potential energy is

\[\begin{equation*}U = mgh = mgl\cos\theta\end{equation*}\]
where \(l\) is the distance from the pivot to the centre of mass. Since we have the components of \(\vec{\omega}\) in the principal axes the rotational kinetic energy about the centre of mass is
\[\begin{align*} T_\mathrm{rot} &= \frac{1}{2} \vec{\omega}\cdot I \vec{\omega} = \frac{1}{2}\left(I_1 \omega_1^2+I_2 \omega_2^2+I_3 \omega_3^2\right) \\ &= \frac{1}{2}\left[ I_1(\dot{\phi}^2\sin^2\theta+\dot{\theta}^2)+ I_3(\dot{\psi}+\dot{\phi}\cos\theta)^2 \right].\end{align*}\]
The instantaneous velocity of the centre of mass can be resolved into two orthogonal directions,
\[\begin{equation*}v_\mathrm{cm} = \begin{pmatrix} l \dot{\theta} \\ l \sin\theta \dot{\phi} \end{pmatrix},\end{equation*}\]
since \(\dot{\vec{\phi}}\) which lies along the \(\hat{z}\) axis is always orthogonal to \(\dot{\vec{\theta}}\) which lies in the \((\hat{x},\hat{y})\) plane. The kinetic energy due to the movement of the centre of mass is therefore
\[\begin{equation}T_\mathrm{cm} = \frac{1}{2} m l^2 \left(\dot{\theta}^2+\sin^2\theta\dot{\phi}^2\right).\end{equation}\]

The total kinetic energy of the top is

\[\begin{equation*}\frac{1}{2}\left[ (I_1+ml^2)(\dot{\phi}^2\sin^2\theta+\dot{\theta}^2)+ I_3(\dot{\psi}+\dot{\phi}\cos\theta)^2 \right].\end{equation*}\]
Notice that the term \(I_1+ml^2=I_1'\) is just the application of the parallel axis theorem to shift the moment of inertia from the centre of mass to the pivot point!

So the Lagrangian is

\[\begin{equation*}\mathcal{L} = \frac{1}{2}\left[ I_1'(\dot{\phi}^2\sin^2\theta+\dot{\theta}^2)+ I_3(\dot{\psi}+\dot{\phi}\cos\theta)^2 \right] -mgl\cos\theta\end{equation*}\]

Now we have the Lagrangian we can derive the equations of motion via the Euler-Lagrange equations. For \(\psi\) we have

\[\begin{equation*}\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot{\psi}} = \frac{\partial \mathcal{L}}{\partial\psi} \;\;\Rightarrow\;\; \frac{d}{dt}\left[I_3(\dot{\psi}+\dot{\phi}\cos\theta)\right] = 0 \;\;\Rightarrow\;\; \frac{d}{dt}\left[I_3\omega_3\right] = 0 \;\;\Rightarrow\;\; \omega_3 = \mathrm{const}\end{equation*}\]
so the `spin’ \(\omega_3\) is conserved. For \(\phi\) we have
\[\begin{gather}\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot{\phi}} = \frac{\partial \mathcal{L}}{\partial\phi} \nonumber \\ \frac{d}{dt}\left(I_3\omega_3\cos\theta +I_1'\dot{\phi}\sin^2\theta\right) = 0 \label{eq:Lzderivative}\\ I_1' \ddot{\phi}\sin\theta -I_3\omega_3\dot{\theta}+2I_1'\dot{\phi}\dot{\theta}\cos\theta = 0\end{gather}\]
where in the last step we have assumed \(\theta\ne 0\) (i.e. top is not perfectly vertical). The \(z\)-component of the angular momentum can be found using \eqref{eq:L} and \eqref{eq:z-axis}:
\[\begin{equation*}L_z = \vec{L}\cdot\hat{{z}} = I_3\omega_3\cos\theta+I_1' \dot{\phi}\sin^2\theta ,\end{equation*}\]
but this is exactly the quantity that is conserved in equation \eqref{eq:Lzderivative}! The last Euler-Lagrange equation is for \(\theta\):
\[\begin{gather}\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot{\theta}} = \frac{\partial \mathcal{L}}{\partial\theta} \nonumber \\ I_1'\ddot{\theta} = \left(I_1'\dot{\phi}^2\cos\theta - I_3\omega_3\dot{\phi} + mgl \right)\sin\theta.\end{gather}\]
It’s usual to label \(\dot{\phi}=\Omega\), which is the rate of precession of the top.

Summarising the results, a symmetrical heavy top is governed by the equations of motion

\[\begin{gather*}I_1' \dot{\Omega}\sin\theta -I_3\omega_3\dot{\theta}+2I_1'\Omega\dot{\theta}\cos\theta = 0 \\ I_1'\ddot{\theta} = \left(I_1'\Omega^2\cos\theta - I_3\omega_3\Omega + mgl \right)\sin\theta,\end{gather*}\]
where \(\omega_3\) is a constant rate of spin and \(\Omega=\dot{\phi}\) is the precession of the top around the \(\hat{{z}}\) axis.

These equations are sufficient to get a solution by numerical integration and this is done below. Notice that the top tends to oscillate in \(\theta\), this is called nutation and we explore it further in the next section.

Click figure to download the CDF demo.
Plot of the intersection of the \(\hat{x}_3\) axis and a unit sphere.

4 Conserved quantities

Three quantities are conserved by this system: the total energy \(E\), and the angular momentum in the \(\hat{x}_3\) and \(\hat{{z}}\) directions. The latter arise since \(\psi\) and \(\phi\) are cyclic variables in the Lagrangian. So we have

\[\begin{gather}L_3 = I_3\omega_3 = I_3(\dot{\psi}+\dot{\phi}\cos\theta) \label{eq:L3conserved},\\ L_z = I_3\omega_3\cos\theta + I_1'\dot{\phi}\sin^2\theta \label{eq:Lzconserved},\\ E = \frac{1}{2}I_3\omega_3^2 +\frac{1}{2}I_1'(\dot{\theta}^2+\dot{\phi}^2\sin^2\theta) + m g l \cos\theta \label{eq:Econserved}.\end{gather}\]
Equations \eqref{eq:L3conserved} and \eqref{eq:Lzconserved} can be rearranged to obtain expressions for \(\dot{\psi}\) and \(\dot{\phi}\) which can then be substituted into \eqref{eq:Econserved} to get,
\[\begin{align}E &= \frac{1}{2}I_1'\dot{\theta}^2 + U(\theta),\label{eq:EandU}\\ U(\theta) &= \frac{1}{2}I_3\omega_3^2 +mgl\cos\theta + \frac{(L_z-I_3\omega_3\cos\theta)^2}{2I_1'\sin^2\theta}.\end{align}\]
This looks just like the equation for the total energy of a particle on a ring subject to the potential \(U(\theta)\). Plotting the potential will give us immediate insight on the kind of solutions that might be available

Click figure to download the CDF demo.
Plot of the nutation potential \(U(\theta)\). The total energy \(E\) is the straight line. Valid \(\theta\) are when \(E\ge U(\theta)\).

From the simulation above it’s clear that in general \(\theta\) is confined to a range \(\theta_\mathrm{min}\le\theta\le\theta_\mathrm{max}\) where \(E\ge U(\theta)\). If \(\omega_3=L_z/I_3\) then \(\theta_\mathrm{min}=0\). This behavious can be readily seen in the numerically integrated solutions in the first simulation.

Finally, we could separate variables and integrate Eq.\eqref{eq:EandU} to get

\[\begin{equation*}t = \int \sqrt{\frac{I_1'}{2[E-U(\theta)]}} d\theta\end{equation*}\]
which is an elliptical integral and at that point we have to resort to numerics anyway.

© Copyright 2022 Alexei Gilchrist