newtonian mechanics - Stability of rotation of a rectangular prism

I've noticed something curious about the rotation of a rectangular prism. If I take a box with height $\neq$ width $\neq$ depth and flip it into the air around different axes of rotation, some motions seem more stable than others. The 3 axes which best illustrate what I mean are:

(1) Through the centre of mass, parallel to the longest box edge.

(2) Through the centre of mass, parallel to the shortest box edge.

(3) Through the centre of mass, parallel to the remaining box edge.

It's "easy" to get the box to rotate cleanly around (1) and (2), but flipping the box around (3) usually results in extra twisting besides the rotation around (3) that I'm trying to achieve (obviously a "perfect" flip on my part would avoid this twisting, which is why I call it an instability). If you're not quite sure what I'm talking about, grab a box or book with 3 different side lengths and try it out (but careful not to break anything!).

What's special about axis (3)?

Image taken from Marsden and Ratiu.

Answer

The rectangular prism is a rigid body. The equations of motion of a rigid body around its center of mass are given by: (Please, see for example: Marsden and Ratiu , (page 6).

$$I_1\dot\Omega_1=(I_2-I_3)\Omega_2\Omega_3$$ $$I_2\dot\Omega_2=(I_3-I_1)\Omega_3\Omega_1$$ $$I_3\dot\Omega_3=(I_1-I_2)\Omega_1\Omega_2$$

Where $\Omega_1,_2,_3$ are the angular velocity components around the body axes and $I_1,_2,_3$ are the corresponding moments of inertia.

Given that the moments of inertia are different, we may assume without loss of generality that: $I_1>I_2>I_3$.

The fact is that the steady motion around the intermediate axis $2$ is not stable, while around the two other axes, the motion is stable. This fact is explained by Marsden and Ratiu on page 30. Also, various other explanations are given in the answers of a related question asked on mathoverflow. Here I'll describe the details of a linearized stability analysis.

A steady state in which the angular velocity vector has only one nonvanishing constant component is a solution of the equations of motion.

For example:

$$\Omega_1=\Omega = const.$$ $$\Omega_2=0$$ $$\Omega_3=0$$

is a solution describing rotation around the first axis. Also

$$\Omega_1=0$$ $$\Omega_2=\Omega = const.$$ $$\Omega_3=0$$

is also a solution describing rotation around the second axis.

Now, we can analyze the stability of small perturbations around these solutions. A perturbation of the first solution is given by:

$$\Omega_1=\Omega + \epsilon \omega_1$$ $$\Omega_2=\epsilon \omega_2$$ $$\Omega_3=\epsilon \omega_3$$

With $\epsilon<<1$. Substituting in the equations of motion and keeping only terms up to the first power of $\epsilon$, we obtain:

$$I_2\dot\omega_2=\epsilon \Omega(I_3-I_1)\omega_3$$ $$I_3\dot\omega_3=\epsilon \Omega(I_1-I_2)\omega_2$$

Taking the first derivative of the second equation with respect to time and substituting the second equation, we obtain:

$$I_2I_3\ddot\omega_3=\epsilon ^2 \Omega^2 (I_3-I_1)(I_1-I_2)\omega_3$$

Since $I_3I_2$, the coefficient on the right hand side is negative and the perturbation satisfies a harmonic oscillator equation of motion of the form:

$$\ddot\omega_3 + k^2 \omega_3 =0$$

Repeating the perturbation analysis for the second solution (rotation about the second axis) we obtain:

$$I_2I_3\ddot\omega_3=\epsilon ^2 \Omega^2 (I_2-I_3)(I_1-I_2)\omega_3$$

Since $I_3I_2$, this coefficient is now negative and the solution describes a Harmonic oscillator with a negative spring constant of the form:

$$\ddot\omega_3 - k^2 \omega_3 =0$$

Which is an unstable perturbation.