2.1   Rotational Kinetic Energy

We consider a rigid body to be made up of a large number of constituent particles with mass m_α, position ${\overrightarrow{x}}_{\alpha }$, and velocities ${\dot{\overrightarrow{x}}}_{\alpha }$, with rigid positional constraints among them. The kinetic energy is

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }\frac{1}{2}{m}_{\alpha }{\dot{\overrightarrow{x}}}_{\alpha }·{\dot{\overrightarrow{x}}}_{\alpha }.}& \left(2.1\right)\end{array}$$

It turns out that the kinetic energy of a rigid body can be separated into two pieces: a kinetic energy of translation and a kinetic energy of rotation. Let’s see how this comes about.

The configuration of a rigid body is fully specified given the location of any point in the body and the orientation of the body. This suggests that it would be useful to decompose the position vectors for the constituent particles as the sum of the vector $\overrightarrow{X}$ to some reference position in the body and the vector ${\overrightarrow{\xi }}_{\alpha }$ from the reference position to the particular constituent element with index α:

$$\begin{array}{cc}{\overrightarrow{x}}_{\alpha }=\overrightarrow{X}+{\overrightarrow{\xi }}_{\alpha }.& \left(2.2\right)\end{array}$$

Along paths, the velocities are related by

$$\begin{array}{cc}{\dot{\overrightarrow{x}}}_{\alpha }=\dot{\overrightarrow{X}}+{\dot{\overrightarrow{\xi }}}_{\alpha }.& \left(2.3\right)\end{array}$$

So in terms of $\dot{\overrightarrow{X}}$ and ${\dot{\overrightarrow{\xi }}}_{\alpha }$ the kinetic energy is

$$\begin{array}{ll}{\displaystyle \sum _{\alpha }\frac{1}{2}{m}_{\alpha }\left(\dot{\overrightarrow{X}}+{\dot{\overrightarrow{\xi }}}_{\alpha }\right)·\left(\dot{\overrightarrow{X}}+{\dot{\overrightarrow{\xi }}}_{\alpha }\right)}\hfill & \hfill \\ ={\displaystyle \sum _{\alpha }\frac{1}{2}{m}_{\alpha }}\left(\dot{\overrightarrow{X}}·\dot{\overrightarrow{X}}+2\dot{\overrightarrow{X}}·{\dot{\overrightarrow{\xi }}}_{\alpha }+{\dot{\overrightarrow{\xi }}}_{\alpha }·{\dot{\overrightarrow{\xi }}}_{\alpha }\right).\hfill & \left(2.4\right)\hfill \end{array}$$

If we select the reference position in the body to be its center of mass,

$$\begin{array}{cc}\overrightarrow{X}=\frac{1}{M}{\displaystyle \sum _{\alpha }{m}_{\alpha }{\overrightarrow{x}}_{\alpha },}& \left(2.5\right)\end{array}$$

where $M={\displaystyle {\sum }_{\alpha }{m}_{\alpha }}$ is the total mass of the body, then

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }{m}_{\alpha }{\overrightarrow{\xi }}_{\alpha }={\displaystyle \sum _{\alpha }{m}_{\alpha }\left({\overrightarrow{x}}_{\alpha }-\overrightarrow{X}\right)=0.}}& \left(2.6\right)\end{array}$$

So along paths the relative velocities satisfy

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }{m}_{\alpha }{\dot{\overrightarrow{\xi }}}_{\alpha }=0.}& \left(2.7\right)\end{array}$$

The kinetic energy is then

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }\frac{1}{2}{m}_{\alpha }\dot{\overrightarrow{X}}·\dot{\overrightarrow{X}}+{\displaystyle {\sum }_{\alpha }\frac{1}{2}{m}_{\alpha }{\dot{\overrightarrow{\xi }}}_{\alpha }·{\dot{\overrightarrow{\xi }}}_{\alpha }.}}& \left(2.8\right)\end{array}$$

The kinetic energy is the sum of the kinetic energy of the motion of the total mass at the center of mass

$$\begin{array}{cc}\frac{1}{2}M\dot{\overrightarrow{X}}·\dot{\overrightarrow{X}},& \left(2.9\right)\end{array}$$

and the kinetic energy of rotation about the center of mass

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }\frac{1}{2}{m}_{\alpha }{\dot{\overrightarrow{\xi }}}_{\alpha }·{\dot{\overrightarrow{\xi }}}_{\alpha }.}& \left(2.10\right)\end{array}$$

Written in terms of appropriate generalized coordinates, the kinetic energy is a Lagrangian for a free rigid body. If we choose generalized coordinates so that the center of mass position is entirely specified by some of them and the orientation is entirely specified by others, then the Lagrange equations for a free rigid body will decouple into two groups of equations, one concerned with the motion of the center of mass and one concerned with the orientation.

Such a separation might occur in other problems, such as a rigid body moving in a uniform gravitational field, but in general, potential energies cannot be separated as the kinetic energy separates. So the motion of the center of mass and the rotational motion are usually coupled through the potential. Even in these cases, it is usually an advantage to choose generalized coordinates that separately specify the position of the center of mass and the orientation.

2.2   Kinematics of Rotation

The motion of a rigid body about a center of rotation, a reference position that is fixed with respect to the body, is characterized at each moment by a rotation axis and a rate of rotation. Let’s elaborate.

We can get from any orientation of a body to any other orientation of the body by a rotation of the body. That this is true is called Euler’s theorem on rotations about a point.² We know that rotations have the property that they do not commute: the composition of successive rotations in general depends on the order of operation. Rotating a book about the $\hat{x}$ axis and then about the ẑ axis puts the book in a different orientation than rotating the book about the ẑ axis and then about the $\hat{x}$ axis. Nevertheless, Euler’s theorem states that however many rotations have been composed to reach a given orientation, the orientation could have been reached with a single rotation. Try it! We take a book, rotate it this way, then that, and then some other way—then find the rotation that does the job in one step. So a rotation can be specified by an axis of rotation and the angular amount of the rotation.

We can specify the orientation of a body by specifying the rotation that takes the body to its actual orientation from some reference orientation. As the body moves, the rotation that does this changes.

Let q be the coordinate path that we will use to describe the motion of the body. Let $\mathcal{M}(q(t))$ be the rotation that takes the body from the reference orientation to the orientation specified by q(t) (see figure 2.1). Let ${\overrightarrow{\xi }}_{\alpha }\left(t\right)$ be the vector to some constituent particle with the body in the orientation specified by q(t), and let $\overrightarrow{{\xi }_{\alpha }^{\prime }}$ be the vector to the same constituent with the body in the reference orientation. Then

$$\begin{array}{cc}{\overrightarrow{\xi }}_{\alpha }\left(t\right)=\mathcal{M}\left(q\left(t\right)\right)\overrightarrow{{\xi }_{\alpha }^{\prime }}.& \left(2.11\right)\end{array}$$

The constituent vectors $\overrightarrow{{\xi }_{\alpha }^{\prime }}$ do not depend on the configuration, because they are the vectors to the positions of the constituents with the body in a fixed reference orientation.

To compute the kinetic energy we will accumulate the contributions from all of the constituent mass elements. So we need the velocities of the constituents. The positions of the constituent particles, at a given time t, are

$$\begin{array}{cc}{\overrightarrow{\xi }}_{\alpha }\left(t\right)=\mathcal{M}\left(q\left(t\right)\right)\overrightarrow{{\xi }_{\alpha }^{\prime }}=M\left(t\right)\overrightarrow{{\xi }_{\alpha }^{\prime }},& \left(2.12\right)\end{array}$$

where $M=\mathcal{M}\circ q$. The velocity is the time derivative

$$\begin{array}{cc}D{\overrightarrow{\xi }}_{\alpha }\left(t\right)=DM\left(t\right)\overrightarrow{{\xi }_{\alpha }^{\prime }}.& \left(2.13\right)\end{array}$$

Using equation (2.12), we can write

$$\begin{array}{cc}D{\overrightarrow{\xi }}_{\alpha }\left(t\right)=DM\left(t\right){\left(M\left(t\right)\right)}^{-1}{\overrightarrow{\xi }}_{\alpha }\left(t\right).& \left(2.14\right)\end{array}$$

So we have a time-varying linear differential equation that describes the motion of the constituents. Let’s look at the multiplier DM(t)(M(t))⁻¹. Since M(t) is a rotation its matrix representation is an orthogonal matrix M(t), with the property ${(\mathbf{M}(t))}^{-1}={(\mathbf{M}(t))}^{\mathcal{T}}$. Because $M(t){(M(t))}^{\mathcal{T}}=I$, its derivative is:

$$\begin{array}{cc}\mathbf{0}=D\left(M{M}^{\mathcal{T}}\right)=DM{M}^{\mathcal{T}}+MD{M}^{\mathcal{T}}.& \left(2.15\right)\end{array}$$

So

$$\begin{array}{cc}DM{M}^{\mathcal{T}}=-{\left(DM{M}^{\mathcal{T}}\right)}^{\mathcal{T}}.& \left(2.16\right)\end{array}$$

We can conclude that $DM{M}^{\mathcal{T}}$ is antisymmetric.

Let u have components (x, y, z). Every 3 × 3 antisymmetric matrix is of the following form:

$$\begin{array}{cc}\mathcal{A}\left(u\right)=\left(\begin{array}{ccc}0& -z& y\\ z& 0& -x\\ -y& x& 0\end{array}\right).& \left(2.17\right)\end{array}$$

Multiplication by this matrix can be interpreted as the operation of cross product with the vector $\overrightarrow{u}$. The vector $\overrightarrow{u}$ has a matrix representation u.

The inverse of the function $\mathcal{A}$ can be applied to any skew-symmetric matrix: we can use ${\mathcal{A}}^{-1}$ to extract the components of u.

We can interpret multiplication by $DM{M}^{\mathcal{T}}$ as a cross product with a vector that we call $\overrightarrow{\omega }$, the angular velocity vector with components ω. So we can write

$$\begin{array}{cc}\mathit{\omega }={\mathcal{A}}^{-1}\left(DM{M}^{\mathcal{T}}\right).& \left(2.18\right)\end{array}$$

In terms of the angular velocity vector, the differential equations for the motion of the constituents (see equation 2.14) are

$$\begin{array}{cc}D{\overrightarrow{\xi }}_{\alpha }\left(t\right)=\overrightarrow{\omega }\left(t\right)\times {\overrightarrow{\xi }}_{\alpha }\left(t\right).& \left(2.19\right)\end{array}$$

If the angular velocity vector for a body is $\overrightarrow{\omega }$ then the velocities of the constituent particles are perpendicular to the vectors to the constituent particles and proportional to the rate of rotation of the body and the distance of the constituent particle from the instantaneous rotation axis:

$$\begin{array}{cc}{\dot{\overrightarrow{\xi }}}_{\alpha }=\overrightarrow{\omega }\times {\overrightarrow{\xi }}_{\alpha }.& \left(2.20\right)\end{array}$$

The components ω′ of the angular velocity vector on the body axes are $\mathit{\omega }\prime =M^{\mathcal{T}}\mathit{\omega }$, so

$$\begin{array}{cc}\mathit{\omega }\prime =M^{\mathcal{T}}{\mathcal{A}}^{-1}\left(DM{M}^{\mathcal{T}}\right).& \left(2.21\right)\end{array}$$

The relationship of the angular velocity vector to the path is a kinematic relationship; it is valid for any path. Thus we can abstract it to obtain the components of the angular velocity at a moment given the configuration and velocity at that moment.

2.3   Moments of Inertia

The rotational kinetic energy is the sum of the kinetic energy of each of the constituents of the rigid body. We can rewrite the rotational kinetic energy in terms of the angular velocity vector and certain aggregate quantities determined by the distribution of mass in the rigid body.

Substituting our representation of the relative velocity vectors into the rotational kinetic energy, we obtain

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }\frac{1}{2}{m}_{\alpha }{\dot{\overrightarrow{\xi }}}_{\alpha }·{\dot{\overrightarrow{\xi }}}_{\alpha }={\displaystyle {\sum }_{\alpha }\frac{1}{2}{m}_{\alpha }\left(\overrightarrow{\omega }\times {\overrightarrow{\xi }}_{\alpha }\right)·\left(\overrightarrow{\omega }\times {\overrightarrow{\xi }}_{\alpha }\right).}}& \left(2.22\right)\end{array}$$

We introduce an arbitrary spatially fixed rectangular coordinate system with origin at the center of rotation and with basis vectors ê₀, ê₁, and ê₂, with the property that ê₀ × ê₁ = ê₂. The components of $\overrightarrow{\omega }$ on this coordinate system are ω⁰, ω¹, and ω². Rewriting $\overrightarrow{\omega }$ in terms of its components, the rotational kinetic energy becomes

$$\begin{array}{c}\begin{array}{l}{\displaystyle \sum _{\alpha }\frac{1}{2}{m}_{\alpha }\left(\left({\displaystyle {\sum }_i{\hat{e}}_i{\omega }^i}\right)\times {\overrightarrow{\xi }}_{\alpha }\right)·\left(\left({\displaystyle {\sum }_j{\hat{e}}_j{\omega }^j}\right)\times {\overrightarrow{\xi }}_{\alpha }\right)}\\ =\frac{1}{2}{\displaystyle {\sum }_{ij}{\omega }^i{\omega }^j}{\displaystyle {\sum }_{\alpha }{m}_{\alpha }}\left({\hat{e}}_i\times {\overrightarrow{\xi }}_{\alpha }\right)·\left({\hat{e}}_j\times {\overrightarrow{\xi }}_{\alpha }\right)\\ =\frac{1}{2}{\displaystyle {\sum }_{ij}{\omega }^i{\omega }^j{I}_{ij},}& \left(2.23\right)\end{array}\end{array}$$

with

$$\begin{array}{cc}{I}_{ij}={\displaystyle \sum _{\alpha }{m}_{\alpha }\left({\hat{e}}_i\times {\overrightarrow{\xi }}_{\alpha }\right)·\left({\hat{e}}_j\times {\overrightarrow{\xi }}_{\alpha }\right).}& \left(2.24\right)\end{array}$$

The nine time-dependent quantities I_ij are the components of the inertia tensor with respect to the chosen coordinate system.

Note what a remarkable form the kinetic energy has taken. All we have done is interchange the order of summations, but now the kinetic energy is written as a sum of products of components of the angular velocity vector, which completely specify how the orientation of the body is changing, and the quantity I_ij, which depends solely on the distribution of mass in the body relative to the chosen coordinate system.

We will deduce a number of properties of the inertia tensor. First, we find a somewhat simpler expression for it. The components of the vector ${\overrightarrow{\xi }}_{\alpha }$ are $\left({\xi }_{\alpha }^0,{\xi }_{\alpha }^1,{\xi }_{\alpha }^2\right)$. If we rewrite ${\overrightarrow{\xi }}_{\alpha }$ as a sum over its components and simplify the elementary vector products of basis vectors, we can obtain the components of the inertia tensor. We can arrange the components of the inertia tensor to form the inertia matrix:

$$\begin{array}{cc}I=\left(\begin{array}{ccc}{I}_{00}& I_{01}& I_{02}\\ I_{10}& I_{11}& I_{12}\\ I_{20}& I_{21}& I_{22}\end{array}\right),& \left(2.25\right)\end{array}$$

where

$$\begin{array}{lll}{I}_{00}\hfill & ={\displaystyle \sum _{\alpha }{m}_{\alpha }({({\xi }_{\alpha }^1)}^2+{({\xi }_{\alpha }^2)}^2)}\hfill & \hfill \\ I_{11}\hfill & ={\displaystyle \sum _{\alpha }{m}_{\alpha }({({\xi }_{\alpha }^0)}^2+{({\xi }_{\alpha }^2)}^2)}\hfill & \hfill \\ I_{22}\hfill & ={\displaystyle \sum _{\alpha }{m}_{\alpha }({({\xi }_{\alpha }^0)}^2+{({\xi }_{\alpha }^1)}^2)}\hfill & \hfill \\ I_{ij}\hfill & =-{\displaystyle \sum _{\alpha }{m}_{\alpha }{\xi }_{\alpha }^i{\xi }_{\alpha }^j}\text{for}i\ne j\hfill & (2.26)\hfill \end{array}$$

Note that the inertia tensor has real components and is symmetric: I_jk = I_kj.

We define the moment of inertia about a line by

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }{m}_{\alpha }{\left({\xi }_{\alpha }^{\perp }\right)}^2},& \left(2.27\right)\end{array}$$

where ${\xi }_{\alpha }^{\perp }$ is the perpendicular distance from the line to the constituent with index α. The diagonal components of the inertia tensor I_ii are recognized as the moments of inertia about the lines coinciding with the coordinate axes ê_i. The off-diagonal components of the inertia tensor are called products of inertia.

The rotational kinetic energy of a body depends on the distribution of mass of the body solely through the inertia tensor. Remarkably, the inertia tensor involves only second-order moments of the mass distribution with respect to the center of mass. We might have expected the kinetic energy to depend in a complicated way on all the moments of the mass distribution, interwoven in some complicated way with the components of the angular velocity vector, but this is not the case. This fact has a remarkable consequence: for the motion of a free rigid body the detailed shape of the body does not matter. If a book and a banana have the same inertia tensor, that is, the same second-order mass moments, then if they are thrown in the same way the subsequent motion will be the same, however complicated that motion is. The facts that the book has corners and the banana has a stem do not affect the motion except for their contributions to the inertia tensor. In general, the potential energy of an extended body is not so simple and does indeed depend on all moments of the mass distribution, but for the kinetic energy the second moments are all that matter!

Exercise 2.1: Rotational kinetic energy

Show that the rotational kinetic energy can also be written

$$\begin{array}{cc}{T}_{\mathrm{R}}=\frac{1}{2}I{\omega }^2,& \left(2.28\right)\end{array}$$

where I is the moment of inertia about the line through the center of mass with direction $\hat{\omega }$, and ω is the instantaneous rate of rotation.

Exercise 2.2: Steiner’s theorem

Let I be the moment of inertia of a body with respect to some given line through the center of mass. Show that the moment of inertia I′ with respect to a second line parallel to the first is

$$\begin{array}{cc}I\prime =I+M{R}^2& \left(2.29\right)\end{array}$$

where M is the total mass of the body and R is the distance between the lines.

Exercise 2.3: Some useful moments of inertia

Show that the moments of inertia of the following objects are as given:

a. The moment of inertia of a sphere of uniform density with mass M and radius R about any line through the center is $\frac{2}{5}M{R}^2$.

b. The moment of inertia of a spherical shell with mass M and radius R about any line through the center is $\frac{2}{3}M{R}^2$.

c. The moment of inertia of a cylinder of uniform density with mass M and radius R about the axis of the cylinder is $\frac{1}{2}M{R}^2$.

d. The moment of inertia of a thin rod of uniform density per unit length with mass M and length L about an axis perpendicular to the rod through the center of mass is $\frac{1}{12}M{L}^2$.

Exercise 2.4: Jupiter

a. The density of a planet increases toward the center. Provide an argument that the moment of inertia of a planet is less than that of a sphere of uniform density of the same mass and radius.

b. The density as a function of radius inside Jupiter is well approximated by

$$\rho \left(r\right)=\frac{M}{R^3}\frac{\sin \left(\pi r/R\right)}{4r/R},$$

where M is the mass and R is the radius of Jupiter. Find the moment of inertia of Jupiter in terms of M and R.

2.4   Inertia Tensor

The representation of the rotational kinetic energy in terms of the inertia tensor was derived with the help of a rectangular coordinate system with basis vectors ê_i. There was nothing special about this particular rectangular basis. So, the kinetic energy must have the same form in any rectangular coordinate system. We can use this fact to derive how the inertia tensor changes if the body or the coordinate system is rotated.

Let’s talk a bit about active and passive rotations. The rotation of the vector $\overrightarrow{x}$ by the rotation R produces a new vector $\overrightarrow{x}\prime =R\overrightarrow{x}$. We may write $\overrightarrow{x}$ in terms of its components with respect to some arbitrary rectangular coordinate system with orthonormal basis vectors ê_i: $\overrightarrow{x}$ = x⁰ê₀ + x¹ê₁ + x²ê₂. Let x indicate the column matrix of components x⁰, x¹, and x² of $\overrightarrow{x}$, and R be the matrix representation of R with respect to the same basis. In these terms the rotation can be written x′ = Rx. The rotation matrix R is a real orthogonal matrix.³ A rotation that carries vectors to new vectors is called an active rotation.

Alternatively, we can rotate the coordinate system by rotating the basis vectors, but leave other vectors that might be represented in terms of them unchanged. If a vector is unchanged but the basis vectors are rotated, then the components of the vector on the rotated basis vectors are not the same as the components on the original basis vectors. Denote the rotated basis vectors by ${\hat{e}}_i^{\prime }=R{\hat{e}}_i$. The component of a vector along a basis vector is the dot product of the vector with the basis vector. So the components of the vector $\overrightarrow{x}$ along the rotated basis ${\hat{e}}_i^{\prime }$ are ${(x\prime )}^i=\overrightarrow{x}·{\hat{e}}_i^{\prime }=\overrightarrow{x}·\left(R{\hat{e}}_i\right)=\left(R^{-1}\overrightarrow{x}\right)·{\hat{e}}_i$.⁴ Thus the components with respect to the rotated basis elements are the same as the components of the rotated vector $R^{-1}\overrightarrow{x}$ with respect to the original basis. In terms of components, if the vector $\overrightarrow{x}$ has components x with respect to the original basis vectors ê_i, then the components x′ of the same vector with respect to the rotated basis vectors ${\hat{e}}_i^{\prime }$ are x′ = R⁻¹x, or equivalently x = Rx′. A rotation that actively rotates the basis vectors, leaving other vectors unchanged, is called a passive rotation. For a passive rotation the components of a fixed vector change as if the vector were actively rotated by the inverse rotation.

With respect to the rectangular basis ê_i the rotational kinetic energy is written

$$\begin{array}{cc}\frac{1}{2}{\displaystyle {\sum }_{ij}{\omega }^i{\omega }^j{I}_{ij}.}& \left(2.30\right)\end{array}$$

In terms of matrix representations, the kinetic energy is

$$\begin{array}{cc}\frac{1}{2}{\mathit{\omega }}^{\mathcal{T}}\mathbf{I}\mathit{\omega },& \left(2.31\right)\end{array}$$

where ω is the column of components representing $\overrightarrow{\omega }$.⁵ If we rotate the coordinate system by the passive rotation R about the center of rotation, the new basis vectors are ${\hat{e}}_i^{\prime }=R{\hat{e}}_i$. The components ω′ of the vector $\overrightarrow{\omega }$ with respect to the rotated coordinate system satisfy

$$\begin{array}{cc}\mathit{\omega }=R\mathit{\omega }\prime ,& \left(2.32\right)\end{array}$$

where R is the matrix representation of R. The kinetic energy is

$$\begin{array}{cc}\frac{1}{2}{(\mathit{\omega }\prime )}^{\mathcal{T}}{R}^{\mathcal{T}}IR\mathit{\omega }\prime .& \left(2.33\right)\end{array}$$

However, if we had started with the basis ${\hat{e}}_i^{\prime }$, we would have written the kinetic energy directly as

$$\begin{array}{cc}\frac{1}{2}{(\mathit{\omega }\prime )}^{\mathcal{T}}I\prime \mathit{\omega }\prime ,& \left(2.34\right)\end{array}$$

where the components are taken with respect to the ${\hat{e}}_i^{\prime }$ basis. Comparing the two expressions, we see that

$$\begin{array}{cc}I\prime =R^{\mathcal{T}}IR.& \left(2.35\right)\end{array}$$

Thus the inertia matrix transforms by a similarity transformation.⁶

2.5   Principal Moments of Inertia

We can use the transformation properties of the inertia tensor (2.35) to show that there are special rectangular coordinate systems for which the inertia tensor I′ is diagonal, that is, $I_{ij}^{\prime }=0$ for i ≠ j. Let’s assume that I′ is diagonal and solve for the rotation matrix R that does the job. Multiplying both sides of (2.35) on the left by R, we have

$$\begin{array}{cc}RI\prime =IR.& \left(2.36\right)\end{array}$$

We can examine pieces of this matrix equation by multiplying on the right by a trivial column vector that picks out a particular column. So we multiply on the right by the column matrix representation e_i of each of the coordinate unit vectors ê_i. These column matrices have a one in the ith row and zeros otherwise. Using $e_i^{\prime }=R{e}_i$, we find

$$\begin{array}{cc}RI\prime e_i=IR{e}_i=I{e}_i^{\prime }.& \left(2.37\right)\end{array}$$

On the other hand, the matrix I′ is diagonal, so

$$\begin{array}{cc}RI\prime e_i=R{e}_i{I}_{ii}^{\prime }=I_{ii}^{\prime }{e}_i^{\prime }.& \left(2.38\right)\end{array}$$

So, from equations (2.37) and (2.38), we have

$$\begin{array}{cc}I{e}_i^{\prime }=I_{ii}^{\prime }{e}_i^{\prime },& \left(2.39\right)\end{array}$$

which we recognize as an equation for the eigenvalue $I_{ii}^{\prime }$ and $e_i^{\prime }$, the column matrix of components of the associated eigenvector.

From $e_i^{\prime }$ = Re_i, we see that the $e_i^{\prime }$ are the columns of the rotation matrix R. Now, rotation matrices are orthogonal, so $R^{\mathcal{T}}R=1$; thus the columns of the rotation matrix must be orthonormal—that is, ${\left(e_i^{\prime }\right)}^{\mathcal{T}}{e}_j^{\prime }={\delta }_{ij}$, where δ_ij is one if i = j and zero otherwise. But the eigenvectors that are solutions of equation (2.39) are not necessarily even orthogonal. So we are not done yet.

If a matrix is real and symmetric then the eigenvalues are real. Furthermore, if the eigenvalues are distinct then the eigenvectors are orthogonal. However, if the eigenvalues are not distinct then the directions of the eigenvectors for the degenerate eigenvalues are not uniquely determined—we have the freedom to choose particular $e_i^{\prime }$ that are orthogonal.⁷ The linearity of equation (2.39) implies that the $e_i^{\prime }$ can be normalized. Thus whether or not the eigenvalues are distinct we can obtain an orthonormal set of $e_i^{\prime }$. This is enough to reconstruct a rotation matrix R that does the job we asked of it: to rotate the coordinate system to a configuration such that the inertia tensor is diagonal. If the eigenvalues are not distinct, the rotation matrix R is not uniquely defined—there is more than one rotation matrix R that does the job.

The eigenvectors and eigenvalues are determined by the requirement that the inertia tensor be diagonal with respect to the rotated coordinate system. Thus the rotated coordinate system has a special orientation with respect to the body. The basis vectors ${\hat{e}}_i^{\prime }$ therefore actually point along particular directions in the body. We define the axes in the body through the center of mass with these directions to be the principal axes. With respect to the coordinate system defined by ${\hat{e}}_i^{\prime }$, the inertia tensor is diagonal, by construction, with the eigenvalues $I_{ii}^{\prime }$ on the diagonal. Thus the moments of inertia about the principal axes are the eigenvalues $I_{ii}^{\prime }$. We call the moments of inertia about the principal axes the principal moments of inertia.

For convenience, we often label the principal moments of inertia according to their size: A ≤ B ≤ C, with principal axis unit vectors â, $\hat{b}$, ĉ, respectively. The positive direction along the principal axes can be chosen so that â, $\hat{b}$, ĉ form a right-handed rectangular coordinate basis.

Let x represent the matrix of components of a vector $\overrightarrow{x}$ with respect to the basis vectors ê_i. Recall that the components x′ of a vector $\overrightarrow{x}$ with respect to the principal axis unit vectors ${\hat{e}}_i^{\prime }$ satisfy

$$\begin{array}{cc}x\prime =R^{\mathcal{T}}x.& \left(2.40\right)\end{array}$$

The components of a vector on the principal axis basis are sometimes called the body components of the vector.

If we choose the reference orientation of the body so that the principal axes are aligned with the spatial axes $\hat{x}$, ŷ, ẑ, then the rotation R that diagonalizes the inertia matrix becomes the rotation M shown in figure 2.1. The axes â′, $\hat{b}\prime$, ĉ′ then become the principal axes. The rotation matrix M multiplies the column of components of a vector on the principal axes to make a column of components of the vector in space.

Now let’s rewrite the kinetic energy in terms of the principal moments of inertia. If we choose our rectangular coordinate system so that it coincides with the principal axes then the calculation is simple. Let the components of the angular velocity vector on the principal axes be (ω^a, ω^b, ω^c). Then, keeping in mind that the inertia tensor is diagonal with respect to the principal axis basis, the kinetic energy is just

$$\begin{array}{cc}{T}_{\mathrm{R}}=\frac{1}{2}\left[A{({\omega }^a)}^2+B{({\omega }^b)}^2+C{({\omega }^c)}^2\right].& (2.41)\end{array}$$

Or as a program:

Exercise 2.5: A constraint on the moments of inertia

Show that the sum of any two of the moments of inertia is greater than or equal to the third moment of inertia. You may assume the moments of inertia are with respect to orthogonal axes.

Exercise 2.6: Principal moments of inertia

For each of the configurations described below find the principal moments of inertia with respect to the center of mass, and find the corresponding principal axes.

a. A regular tetrahedron consisting of four equal point masses tied together with rigid massless wire.

b. A cube of uniform density.

c. Five equal point masses rigidly connected by massless stuff. The point masses are at the rectangular coordinates

(−1, 0, 0), (1, 0, 0), (1, 1, 0), (0, 0, 0), (0, 0, 1).

Exercise 2.7: This book

Measure this book. You will admit that it is pretty dense. Don’t worry, you will get to throw it later. Show that the principal axes are the lines connecting the centers of opposite faces of the idealized brick approximating the book. Compute the corresponding principal moments of inertia.

2.6   Vector Angular Momentum

The vector angular momentum of a particle is the cross product of its position vector and its linear momentum vector. For a rigid body the vector angular momentum is the sum of the vector angular momentum of each of the constituents. Here we find an expression for the vector angular momentum of a rigid body in terms of the inertia tensor and the angular velocity vector.

The vector angular momentum of a rigid body is

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }{\overrightarrow{x}}_{\alpha }\times \left(m_{\alpha }{\dot{\overrightarrow{x}}}_{\alpha }\right),}& \left(2.42\right)\end{array}$$

where ${\overrightarrow{x}}_{\alpha }$, ${\dot{\overrightarrow{x}}}_{\alpha }$, and m_α are the positions, velocities, and masses of the constituent particles. It turns out that the vector angular momentum decomposes into the sum of the angular momentum of the center of mass and the rotational angular momentum about the center of mass, just as the kinetic energy separates into the kinetic energy of the center of mass and the kinetic energy of rotation. As in the kinetic energy demonstration (section 2.1), decompose the position into the vector to the center of mass $\overrightarrow{X}$ and the vectors from the center of mass to the constituent mass elements ${\overrightarrow{\xi }}_{\alpha }$:

$$\begin{array}{cc}{\overrightarrow{x}}_{\alpha }=\overrightarrow{X}+{\overrightarrow{\xi }}_{\alpha },& \left(2.43\right)\end{array}$$

with velocities

$$\begin{array}{cc}{\dot{\overrightarrow{x}}}_{\alpha }=\dot{\overrightarrow{X}}+{\dot{\overrightarrow{\xi }}}_{\alpha }.& \left(2.44\right)\end{array}$$

Substituting, the angular momentum is

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }{m}_{\alpha }\left(\overrightarrow{X}+{\overrightarrow{\xi }}_{\alpha }\right)\times \left(\dot{\overrightarrow{X}}+{\dot{\overrightarrow{\xi }}}_{\alpha }\right).}& \left(2.45\right)\end{array}$$

Multiplying out the product, and using the fact that $\overrightarrow{X}$ is the center of mass and $M={\displaystyle {\sum }_{\alpha }{m}_{\alpha }}$ is the total mass of the body, the angular momentum is

$$\begin{array}{cc}\overrightarrow{X}\times \left(M\dot{\overrightarrow{X}}\right)+{\displaystyle \sum _{\alpha }{\overrightarrow{\xi }}_{\alpha }\times \left(m_{\alpha }{\dot{\overrightarrow{\xi }}}_{\alpha }\right).}& \left(2.46\right)\end{array}$$

The angular momentum of the center of mass is

$$\begin{array}{cc}\overrightarrow{X}\times \left(M\dot{\overrightarrow{X}}\right),& \left(2.47\right)\end{array}$$

and the rotational angular momentum is

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }{\overrightarrow{\xi }}_{\alpha }\times \left(m_{\alpha }{\dot{\overrightarrow{\xi }}}_{\alpha }\right).}& \left(2.48\right)\end{array}$$

Using ${\dot{\overrightarrow{\xi }}}_{\alpha }=\overrightarrow{\omega }\times {\overrightarrow{\xi }}_{\alpha }$, we get the rotational angular momentum vector

$$\begin{array}{cc}\overrightarrow{L}={\displaystyle \sum _{\alpha }{m}_{\alpha }{\overrightarrow{\xi }}_{\alpha }\times \left(\overrightarrow{\omega }\times {\overrightarrow{\xi }}_{\alpha }\right).}& \left(2.49\right)\end{array}$$

We can also reexpress the rotational angular momentum in terms of the angular velocity vector and the inertia tensor, as we did for the kinetic energy. In terms of components with respect to the basis {ê₀, ê₁, ê₂}, this is

$$\begin{array}{cc}{L}_j={\displaystyle \sum _k{I}_{jk}{\omega }^k},& \left(2.50\right)\end{array}$$

where I_jk are the components of the inertia tensor (2.24). The angular momentum and the kinetic energy are expressed in terms of the same inertia tensor.

With respect to the principal-axis basis, the components of the angular momentum have a particularly simple form:

$$\begin{array}{lll}{L}_a\hfill & =A{\omega }^a\hfill & \hfill \\ L_b\hfill & =B{\omega }^b\hfill & \hfill \\ L_c\hfill & =C{\omega }^c\hfill & \left(2.51\right)\hfill \end{array}$$

Since the angular momenta are the partial derivatives of T_R (see equation 2.41) with respect to the angular velocities, they must be grouped as a down tuple (in matrix language, a row matrix): L′ = [L_a, L_b, L_c]. As a program:

If M is the matrix representation of the rotation that takes an angular-velocity vector $\overrightarrow{\omega }\prime$ to a rotated vector $\overrightarrow{\omega }$, the components transform as ω = Mω′.

When working with matrices it is more convenient to work with a column matrix of the angular momentum components, so we introduce $\overline{L}=L^{\mathcal{T}}$. Using ω = Mω′ and equation (2.35) with R replaced by M we derive an expression for the angular momentum

$$\begin{array}{cc}\overline{L}=I\mathit{\omega }=MI\prime \mathit{\omega }\prime =M\overline{L}\prime .& \left(2.52\right)\end{array}$$

Transposing this result, we see that the angular momentum components must transform as $L=L\prime M^{\mathcal{T}}$:

Exercise 2.8: Rotational angular momentum

Verify that expression (2.50) for the components of the rotational angular momentum (2.49) in terms of the inertia tensor is correct.

2.7   Euler Angles

To go further we must finally specify a set of generalized coordinates. We first do this using the traditional Euler angles. Later, we find other ways of describing the orientation of a rigid body.

We are using an intermediate representation of the orientation in terms of the function $\mathcal{M}$ of the generalized coordinates that gives the rotation that takes the body from some reference orientation and rotates it to the orientation specified by the generalized coordinates. Here we take the reference orientation so that principal-axis unit vectors â, $\hat{b}$, ĉ are coincident with the basis vectors ê_i, labeled here by $\hat{x}$, ŷ, ẑ.

We define the Euler angles in terms of simple rotations about the coordinate axes. Let R_x(ψ) be a right-handed rotation about the $\hat{x}$ axis by the angle ψ, and let R_z(ψ) be a right-handed rotation about the ẑ axis by the angle ψ. The function $\mathcal{M}$ for Euler angles is written as a composition of three of these simple coordinate axis rotations:

$$\begin{array}{cc}\mathcal{M}\left(\theta ,\phi ,\psi \right)=R_z\left(\phi \right)\circ R_x\left(\theta \right)\circ R_z\left(\psi \right),& \left(2.53\right)\end{array}$$

for the Euler angles θ, φ, ψ.

The Euler angles can specify any orientation of the body, but the orientation does not always correspond to a unique set of Euler angles. In particular, if θ = 0 then the orientation is dependent only on the sum φ + ψ, so the orientation does not uniquely determine either φ or ψ.

Exercise 2.9: Euler angles

It is not immediately obvious that all orientations can be represented in terms of the Euler angles. To show that the Euler angles are adequate to represent all orientations, solve for the Euler angles that give an arbitrary rotation R. Keep in mind that some orientations do not correspond to a unique representation in terms of Euler angles.

Though the Euler angles allow us to specify all orientations and thus can be used as generalized coordinates, the definition of Euler angles is pretty arbitrary. In fact no reasoning has led us to them. This is reflected in our presentation of them by just saying “here they are.” Euler angles are well suited for some problems, but cumbersome for others.

There are other ways of defining similar sets of angles. For instance, we could also take our generalized coordinates to satisfy

$$\begin{array}{cc}\mathcal{M}\prime \left(\theta ,\phi ,\psi \right)=R_x\left(\phi \right)\circ R_y\left(\theta \right)\circ R_z\left(\psi \right).& \left(2.54\right)\end{array}$$

Such alternatives to the Euler angles sometimes come in handy.

Each of the fundamental rotations can be represented as a matrix. The rotation matrix representing a right-handed rotation about the ẑ axis by the angle ψ is

$$\begin{array}{cc}{R}_z\left(\psi \right)=\left(\begin{array}{ccc}\cos \psi & -\sin \psi & 0\\ \sin \psi & \cos \psi & 0\\ 0& 0& 1\end{array}\right)& \left(2.55\right)\end{array}$$

and a right-handed rotation about the x axis by the angle ψ is represented by the matrix

$$\begin{array}{cc}{R}_x\left(\psi \right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \psi & -\sin \psi \\ 0& \sin \psi & \cos \psi \end{array}\right).& \left(2.56\right)\end{array}$$

The matrix that represents the rotation that carries the body from its reference orientation to the actual orientation is

$$\begin{array}{cc}M\left(\theta ,\phi ,\psi \right)=R_z\left(\phi \right)R_x\left(\theta \right)R_z\left(\psi \right).& \left(2.57\right)\end{array}$$

The rotation matrices and their product can be constructed by simple programs:

Now that we have a procedure that implements a sample $\mathcal{M}$, we can find the components of the angular velocity vector and the body components of the angular velocity vector using the procedures M-of-q->omega-of-t and M-of-q->omega-body-of-t from section 2.2. For example,

To construct the kinetic energy we need the procedure of state that gives the body components of the angular velocity vector:

We capture this result as a procedure:

The kinetic energy can now be written:

We can define procedures to calculate the components of the angular momentum on the principal axes:

We then transform the components of the angular momentum on the principal axes to the components on the fixed basis ê_i:

These procedures are local state functions, like Lagrangians.

2.8   Motion of a Free Rigid Body

The kinetic energy, expressed in terms of a suitable set of generalized coordinates, is a Lagrangian for a free rigid body. In section 2.1 we found that the kinetic energy of a rigid body can be written as the sum of the rotational kinetic energy and the translational kinetic energy. If we choose one set of coordinates to specify the position and another set to specify the orientation, the Lagrangian becomes a sum of a translational Lagrangian and a rotational Lagrangian. The Lagrange equations for translational motion are not coupled to the Lagrange equations for the rotational motion. For a free rigid body the translational motion is just that of a free particle: uniform motion. Here we concentrate on the rotational motion of the free rigid body. We can adopt the Euler angles as the coordinates that specify the orientation; the rotational kinetic energy was expressed in terms of Euler angles in the previous section.

Conserved quantities

The Lagrangian for a free rigid body has no explicit time dependence, so we can deduce that the energy, which is just the kinetic energy, is conserved by the motion.

The Lagrangian does not depend on the Euler angle φ, so we can deduce that the momentum conjugate to this coordinate is conserved. An explicit expression for the momentum conjugate to φ is

(define Euler-state

(up ’t

(up ’theta ’phi ’psi)

(up ’thetadot ’phidot ’psidot)))

(show-expression

(ref (((partial 2) (T-body-Euler ’A ’B ’C)) Euler-state)

1))

$\begin{array}{l}A\dot{\phi }{\left(\sin \left(\theta \right)\right)}^2{\left(\sin \left(\psi \right)\right)}^2+A\dot{\theta }\cos \left(\psi \right)\sin \left(\theta \right)\sin \left(\psi \right)\hfill \\ +B\dot{\phi }{\left(\cos \left(\psi \right)\right)}^2{\left(\sin \left(\theta \right)\right)}^2-B\dot{\theta }\cos \left(\psi \right)\sin \left(\theta \right)\sin \left(\psi \right)\hfill \\ +C\dot{\phi }{\left(\cos \left(\theta \right)\right)}^2+C\dot{\psi }\cos \left(\theta \right)\hfill \end{array}$

We know that this complicated quantity is conserved by the motion of the rigid body because of the symmetries of the Lagrangian.

If there are no external torques, then we expect that the vector angular momentum will be conserved. We can verify this using the Lagrangian formulation of the problem. First, we note that L_z is the same as p_φ. We can check this by direct calculation:

(- (ref ((L-space-Euler ’A ’B ’C) Euler-state)

2)

(ref (((partial 2) (T-body-Euler ’A ’B ’C)) Euler-state)

1))

0

We know that p_φ is conserved because the Lagrangian for the free rigid body did not mention φ, so now we know that L_z is conserved. Since the orientation of the coordinate axes is arbitrary, we know that if any rectangular component is conserved then all of them are. So the vector angular momentum is conserved for the free rigid body.

We could have seen this with the help of Noether’s theorem (see section 1.8.5). There is a continuous family of rotations that can transform any orientation into any other orientation. The orientation of the coordinate axes we used to define the Euler angles is arbitrary, and the kinetic energy (the Lagrangian) is the same for any choice of coordinate system. Thus the situation meets the requirements of Noether’s theorem, which tells us that there is a conserved quantity. In particular, the family of rotations around each coordinate axis gives us conservation of the angular momentum component on that axis. We construct the vector angular momentum by combining these contributions.

Exercise 2.10: Uniformly accelerated rigid body

Show that a rigid body subject to a uniform acceleration rotates as a free rigid body, while the center of mass has a parabolic trajectory.

Exercise 2.11: Conservation of angular momentum

Fill in the details of the argument that Noether’s theorem implies that vector angular momentum is conserved by the motion of the free rigid body.

2.8.1 Computing the Motion of Free Rigid Bodies

Lagrange’s equations for the motion of a free rigid body in terms of Euler angles are quite disgusting, so we will not show them here. However, we will use the Lagrange equations to explore the motion of the free rigid body.

Before doing this it is worth noting that the equations of motion in Euler angles are singular for some configurations, because for these configurations the Euler angles are not uniquely defined. If we set θ = 0 then an orientation does not correspond to a unique value of φ and ψ; only their sum determines the orientation.

The singularity arises in the explicit Lagrange equations when we attempt to solve for the second derivative of the generalized coordinates in terms of the generalized coordinates and the generalized velocities (see section 1.7). The isolation of the second derivative requires multiplying by the inverse of ∂₂∂₂L. The determinant of this quantity becomes zero when the Euler angle θ is zero:

(show-expression

(determinant

(((square (partial 2)) (T-body-Euler ’A ’B ’C))

Euler-state)))

$$ABC{\left(\sin \left(\theta \right)\right)}^2$$

So when θ is zero, we cannot solve for the second derivatives. When θ is small, the Euler angles can move very rapidly, and thus may be difficult to compute reliably. Of course, the motion of the rigid body is perfectly well behaved for any orientation. This is a problem of the representation of that motion in Euler angles; it is a “coordinate singularity.”

One solution to this problem is to use another set of Euler-like coordinates for which Lagrange’s equations have singularities for different orientations, such as those defined in equation (2.54). So if as the calculation proceeds the trajectory comes close to a singularity in one set of coordinates, we can switch coordinate systems and use another set for a while until the trajectory encounters another singularity. This solves the problem, but it is cumbersome. For the moment we will ignore this problem and compute some trajectories, being careful to limit our attention to trajectories that avoid the singularities.

We will compute some trajectories by numerical integration and check our integration process by seeing how well energy and angular momentum are conserved. Then, we will investigate the evolution of the components of angular momentum on the principal axis basis. We will discover that we can learn quite a bit about the qualitative behavior of rigid bodies by combining the information we get from the energy and angular momentum.

To develop a trajectory from initial conditions we integrate the Lagrange equations, as we did in chapter 1. The system derivative is obtained from the Lagrangian:

(define (rigid-sysder A B C)

(Lagrangian->state-derivative (T-body-Euler A B C)))

The following program monitors the errors in the energy and in the components of the angular momentum:

(define ((monitor-errors win A B C L0 E0) state)

(let ((t (time state))

(L ((L-space-Euler A B C) state))

(E ((T-body-Euler A B C) state)))

(plot-point win t (relative-error (ref L 0) (ref L0 0)))

(plot-point win t (relative-error (ref L 1) (ref L0 1)))

(plot-point win t (relative-error (ref L 2) (ref L0 2)))

(plot-point win t (relative-error E E0))))

(define (relative-error value reference-value)

(if (zero? reference-value)

(error “Zero reference value -- RELATIVE-ERROR”)

(/ (- value reference-value) reference-value)))

We make a plot window to display the errors:

(define win (frame 0.0 100.0 -1.0e-12 1.0e-12))

The default integration method used by the system is Bulirsch–Stoer (bulirsch-stoer), but here we set the integration method to be quality-controlled Runge–Kutta (qcrk4), because the error plot is more interesting:

(set-ode-integration-method! ’qcrk4)

We use evolve to investigate the evolution:

(let ((A 1.0) (B (sqrt 2.0)) (C 2.0)	; moments of inertia
(state0 (up 0.0	; initial state

(up 1.0 0.0 0.0)

(up 0.1 0.1 0.1))))

(let ((L0 ((L-space-Euler A B C) state0))

(E0 ((T-body-Euler A B C) state0)))

((evolve rigid-sysder A B C)

state0

(monitor-errors win A B C L0 E0)

0.1	; step between plotted points
100.0	; final time
1.0e-12)))	; max local truncation error

The plot that is developed of the relative errors in the components of the angular momenta and the energy (see figure 2.2) shows that we have been successful in controlling the error in the conserved quantities. This should give us some confidence in the trajectory that is evolved.

2.8.2 Qualitative Features of Free Rigid Body Motion

The evolution of the components of the angular momentum on the principal axes has a remarkable property. For almost every initial condition the body components of the angular momentum periodically trace a simple closed curve.

We can see this by investigating a number of trajectories and plotting the components of angular momentum of the body on the principal axes (see figure 2.3). To make this figure a number of trajectories of equal energy were computed. The three-dimensional space of body components is projected onto a two-dimensional plane for display. Points on the back of this projection of the ellipsoid of constant energy are plotted with lower density than points on the front of the ellipsoid. For most initial conditions we find a one-dimensional simple closed curve. Some trajectories on the front side appear to cross trajectories on the back side, but this is an artifact of projection. There is also a family of trajectories that appear to intersect in two points, one on the front side and one on the back side. The curve that is the union of these trajectories is called a separatrix; it separates different types of motion.

What is going on? The state space for a free rigid body is six-dimensional: the three Euler angles and their time derivatives. We know four constants of the motion—the three spatial components of the angular momentum, L_x, L_y, and L_z, and the energy, E. Thus, the motion is restricted to a two-dimensional region of the state space.⁸ Our experiment shows that the components of the angular momentum trace one-dimensional closed curves in the angular-momentum subspace, so there is something more going on here.

The total angular momentum is conserved if all of the components are, so we also have the constant

$$\begin{array}{cc}{L}^2=L_x^2+L_y^2+L_z^2.& \left(2.58\right)\end{array}$$

The spatial components of the angular momentum do not change, but of course the projections of the angular momentum onto the principal axes do change because the axes move as the body moves. However, the magnitude of the angular momentum vector is the same whether it is computed from components on the fixed basis or components on the principal axis basis. So, the combination

$$\begin{array}{cc}{L}^2=L_a^2+L_b^2+L_c^2,& \left(2.59\right)\end{array}$$

is conserved.

Using the expressions (2.51) for the components of the angular momentum in terms of the components of the angular velocity vector on the principal axes, the kinetic energy (2.41) can be rewritten in terms of the angular momentum components on the principal axes:

$$\begin{array}{cc}E=\frac{1}{2}\left(\frac{L_a^2}{A}+\frac{L_b^2}{B}+\frac{L_c^2}{C}\right).& \left(2.60\right)\end{array}$$

The two conserved quantities (2.59 and 2.60) provide constraints on how the components of the angular momentum vector on the principal axes can change. We recognize the conservation of angular momentum constraint (2.59) as the equation of a sphere, and the conservation of kinetic energy constraint (2.60) as the equation for a triaxial ellipsoid. For every trajectory both constraints are satisfied, so the components of the angular momentum move on the intersection of these two surfaces, the energy ellipsoid and the angular momentum sphere. The intersection of an ellipsoid and a sphere with the same center is typically two closed curves, so an orbit is confined to one of these curves. This sheds light on the puzzle posed at the beginning of this section.

Because of our ordering A ≤ B ≤ C, the longest axis of this triaxial ellipsoid coincides with the ĉ direction (all the angular momentum is along the axis of largest principal moment of inertia) and the shortest axis of the energy ellipsoid coincides with the â axis (all the angular momentum is along the smallest moment of inertia). Without actually solving the Lagrange equations, we have found strong constraints on the evolution of the components of the angular momentum on the principal axes.

To determine how the system evolves along these intersection curves we have to use the equations of motion. We observe that the evolution of the components of the angular momentum on the principal axes depends only on the components of the angular momentum on the principal axes, even though the values of these components are not enough to completely specify the dynamical state. Apparently the dynamics of these components is self-contained, and we will see that it can be described in terms of a set of differential equations whose only dynamical variables are the components of the angular momentum on the principal axes (see section 2.9).

We note that there are two axes for which the intersection curves shrink to a point if we hold the energy constant and vary the magnitude of the angular momentum. If the angular momentum starts at these points, the conserved quantities constrain the angular momentum to stay there. These points are equilibrium points for the body components of the angular momentum. However, they are not equilibrium points for the system as a whole. At these points the body is still rotating even though the body components of the angular momentum are not changing. This kind of equilibrium is called a relative equilibrium. We can also see that if the angular momentum is initially slightly displaced from one of these relative equilibria, then the angular momentum is constrained to stay near it on one of the intersection curves. The angular momentum vector is fixed in space, so the principal axis of the equilibrium point of the body rotates stably about the angular momentum vector.

At the principal axis with intermediate moment of inertia, the $\hat{b}$ axis, the intersection curves appear to cross. As we observed, the dynamics of the components of the angular momentum on the principal axes forms a self-contained dynamical system. Trajectories of a dynamical system cannot cross,⁹ so the most that can happen is that if the equations of motion carry the system along the intersection curve then the system can approach the crossing point only asymptotically. So without solving any equations we can deduce that the point of crossing is another relative equilibrium. If the angular momentum is initially aligned with the intermediate axis, then it stays aligned. If the system is slightly displaced from the intermediate axis, then the evolution along the intersection curve will take the system far from the relative equilibrium. So rotation about the axis of intermediate moment of inertia is unstable—initial displacements of the angular momentum, however small initially, become large. Again, the angular momentum vector is fixed in space, but now the principal axis with the intermediate principal moment does not stay close to the angular momentum, so the body executes a complicated tumbling motion.

This gives some insight into the mystery of the thrown book mentioned at the beginning of the chapter. If one throws a book so that it is initially rotating about either the axis with the largest moment of inertia or the axis with the smallest moment of inertia (the shortest and longest physical axes, respectively), the book rotates regularly about that axis. However, if the book is thrown so that it is initially rotating about the axis of intermediate moment of inertia (the intermediate physical axis), then it tumbles, however carefully it is thrown. You can try it with this book (but put a rubber band or string around it first).

Before moving on, we can make some further physical deductions. Suppose a freely rotating body is subject to some sort of internal friction that dissipates energy but conserves the angular momentum. For example, real bodies flex as they spin. If the spin axis moves with respect to the body then the flexing changes with time, and this changing distortion converts kinetic energy of rotation into heat. Internal processes do not change the total angular momentum of the system. If we hold the magnitude of the angular momentum fixed but gradually decrease the energy, then the curve of intersection on which the system moves gradually deforms. For a given angular momentum there is a lower limit on the energy: the energy cannot be so low that there are no intersections. For this lowest energy the intersection of the angular momentum sphere and the energy ellipsoid is a pair of points on the axis of maximum moment of inertia. With energy dissipation, a freely rotating physical body eventually ends up with the lowest energy consistent with the given angular momentum, which is rotation about the principal axis with the largest moment of inertia (typically the shortest physical axis).

Thus, we expect that given enough time all freely rotating physical bodies will end up rotating about the axis of largest moment of inertia. You can demonstrate this to your satisfaction by twirling a small bottle containing some viscous fluid, such as correction fluid. What you will find is that, whatever spin you try to put on the bottle, it will reorient itself so that the axis of the largest moment of inertia is aligned with the spin axis. Remarkably, this is very nearly true of almost every body in the solar system for which there is enough information to decide. The deviations from principal axis rotation for the Earth are tiny: the angle between the angular momentum vector and the ĉ axis for the Earth is less than one arc-second.¹⁰ In fact, the evidence is that all of the planets, the Moon and all the other natural satellites, and almost all of the asteroids rotate very nearly about the largest moment of inertia. We have deduced that this is to be expected using an elementary argument. There are exceptions. Comets typically do not rotate about the largest moment. As they are heated by the sun, material spews out from localized jets, and the back reaction from these jets changes the rotation state. Among the natural satellites, the only known exception is Saturn’s satellite Hyperion, which is tumbling chaotically. Hyperion is especially out of round and subject to strong gravitational torques from Saturn.

2.9   Euler’s Equations

For a free rigid body we have seen that the components of the angular momentum on the principal axes comprise a self-contained dynamical system: the variation of the principal axis components depends only on the principal axis components. Here we derive equations that govern the evolution of these components.

The starting point for the derivation is the conservation of the vector angular momentum. The components of the angular momentum on the principal axes are¹¹

$$\begin{array}{cc}\overline{L}\prime =I\prime \mathit{\omega }\prime ,& \left(2.61\right)\end{array}$$

where ω′ is composed of the components of the angular velocity vector on the principal axes and I′ is the matrix representation of the inertia tensor with respect to the principal axis basis:

$$\begin{array}{cc}I\prime =\left(\begin{array}{ccc}A& 0& 0\\ 0& B& 0\\ 0& 0& C\end{array}\right).& \left(2.62\right)\end{array}$$

The body components of the angular momentum L′ are related to the components L on the fixed rectangular basis ê_i by

$$\begin{array}{cc}\overline{L}=\mathbf{M}\overline{L}\prime ,& \left(2.63\right)\end{array}$$

where M is the matrix representation of the rotation that carries the body and all vectors attached to the body from the reference orientation of the body to the actual orientation.

The vector angular momentum is conserved for free rigid-body motion, and so are its components on a fixed rectangular basis. So, along solution paths,

$$\begin{array}{cc}0=D\overline{L}=DM\overline{L}\prime +MD\overline{L}\prime .& \left(2.64\right)\end{array}$$

Solving, we find

$$\begin{array}{cc}D\overline{L}\prime =-M^{\mathcal{T}}DM\overline{L}\prime .& \left(2.65\right)\end{array}$$

In terms of ω′ this is

$$\begin{array}{lll}I\prime D\mathit{\omega }\prime \hfill & =-M^{\mathcal{T}}DMI\prime \mathit{\omega }\prime \hfill & \hfill \\ \hfill & =-M^{\mathcal{T}}\mathcal{A}\left(M\mathit{\omega }\prime \right)MI\prime \mathit{\omega }\prime ,\hfill & \left(2.66\right)\hfill \end{array}$$

where we have used equation (2.21) to write DM in terms of $\mathcal{A}$. The function $\mathcal{A}$ has the property¹²

$$\begin{array}{cc}{R}^{\mathcal{T}}\mathcal{A}\left(Rv\right)R=\mathcal{A}\left(v\right)& \left(2.67\right)\end{array}$$

for any vector with components v and any rotation with matrix representation R. Using this property of $\mathcal{A}$, we find Euler’s equations:

$$\begin{array}{cc}I\prime D\mathit{\omega }\prime =-\mathcal{A}(\mathit{\omega }\prime )I\prime \mathit{\omega }\prime .& \left(2.68\right)\end{array}$$

Euler’s equations give the time derivative of the body components of the angular velocity vector entirely in terms of the angular velocity components and the principal moments of inertia. Let ω^a, ω^b, and ω^c denote the components of the angular velocity vector on the principal axes. Then Euler’s equations can be written as the component equations

$$\begin{array}{ll}AD{\omega }^a=\left(B-C\right){\omega }^b{w}^c\hfill & \hfill \\ BD{\omega }^b=\left(C-A\right){\omega }^c{\omega }^a\hfill & \hfill \\ CD{\omega }^c=\left(A-B\right){\omega }^a{\omega }^b.\hfill & \left(2.69\right)\hfill \end{array}$$

Alternatively, we can rewrite Euler’s equations in terms of the components of the angular momentum on the principal axes

$$\begin{array}{ll}D{L}_a=\left(\frac{1}{C}-\frac{1}{B}\right)L_b{L}_c\hfill & \hfill \\ D{L}_b=\left(\frac{1}{A}-\frac{1}{C}\right)L_a{L}_c\hfill & \hfill \\ D{L}_a=\left(\frac{1}{B}-\frac{1}{A}\right)L_a{L}_b.\hfill & \left(2.70\right)\hfill \end{array}$$

These equations confirm that the time derivatives of the components of the angular momentum on the principal axes depend only on the components of the angular momentum on the principal axes.

Euler’s equations are very simple, but they do not completely determine the evolution of a rigid body—they do not give the spatial orientation of the body. However, equation (2.21) and property (2.67) can be used to relate the derivative of the orientation matrix to the body components of the angular velocity vector:

$$\begin{array}{cc}DM=M\mathcal{A}(\mathit{\omega }\prime ).& \left(2.71\right)\end{array}$$

A straightforward method of using these equations is to integrate them componentwise as a set of nine first-order ordinary differential equations, with initial conditions determining the initial configuration matrix. Together with Euler’s equations, which describe how the body components of the angular velocity vector change with time, this system of equations governing the motion of a rigid body is complete. However, the reader will no doubt have noticed that this approach is rather wasteful. The fact that the orientation matrix can be specified with only three parameters has not been taken into account. We should be integrating three equations for the orientation, given ω′, not nine. To accomplish this we once again need to parameterize the configuration matrix.

For example, we can use Euler angles to parameterize the orientation:

$$\begin{array}{cc}\mathcal{M}\left(\theta ,\phi ,\psi \right)=R_z\left(\phi \right)R_x\left(\theta \right)R_z\left(\psi \right).& \left(2.72\right)\end{array}$$

We form M by composing $\mathcal{M}$ with an Euler coordinate path. Equation (2.71) can then be used to solve for Dθ, Dφ, and Dψ. We find

$$\begin{array}{cc}\left(\begin{array}{l}D\theta \hfill \\ D\phi \hfill \\ D\psi \hfill \end{array}\right)=\frac{1}{\sin \theta }\left(\begin{array}{ccc}\cos \psi \sin \theta & -\sin \psi \sin \theta & 0\\ \sin \psi & \cos \psi & 0\\ -\sin \psi \cos \theta & -\cos \psi \cos \theta & \sin \theta \end{array}\right)\left(\begin{array}{c}{\omega }^a\\ {\omega }^b\\ {\omega }^c\end{array}\right).& \left(2.73\right)\end{array}$$

This gives us the desired equation for the orientation. Note that it is singular for θ = 0, as are Lagrange’s equations. So Euler’s equations using Euler angles for the configuration have the same problem as did the Lagrange equations using Euler angles. Again, this is a manifestation of the fact that for θ = 0 the orientation depends only on φ+ψ. The singularity in the equations of motion for θ = 0 does not correspond to anything funny in the motion of the rigid body. A practical solution to the singularity problem is to choose another set of Euler-like angles that have a singularity in a different place, and switch from one to the other when the going gets tough.

Exercise 2.12:

Fill in the details of the derivation of equation (2.73). You may want to use the computer to help with the algebra.

2.10   Axisymmetric Tops

We have all played with a top at one time or another. For the purposes of analysis we will consider an idealized top that does not wander around. Thus, an ideal top is a rotating rigid body, one point of which is fixed in space. Furthermore, the center of mass of the top is not at the fixed point, which is the center of rotation, and there is a uniform gravitational acceleration.

For our top we can take the Lagrangian to be the difference of the kinetic energy and the potential energy. We already know how to write the kinetic energy—what is new here is that we must express the potential energy in terms of the configuration. In the case of a body in a uniform gravitational field this is easy. The potential energy is the sum of “mgh” for all the constituent particles:

$$\begin{array}{cc}{\displaystyle \sum _{\alpha }{m}_{\alpha }g{h}_{\alpha },}& \left(2.84\right)\end{array}$$

where g is the gravitational acceleration, h_α = ${\overrightarrow{x}}_{\alpha }$ · ẑ, and the unit vector ẑ indicates which way is up. Rewriting the vector to the constituents in terms of the vector $\overrightarrow{X}$ to the center of mass, the potential energy is

$$\begin{array}{ll}{\displaystyle \sum _{\alpha }{m}_{\alpha }g\left(\overrightarrow{X}+{\overrightarrow{\xi }}_{\alpha }\right)·\hat{z}}\hfill & \hfill \\ =gM\overrightarrow{X}·\hat{z}+g\left({\displaystyle \sum _{\alpha }{m}_{\alpha }{\overrightarrow{\xi }}_{\alpha }}\right)·\hat{z}\hfill & \hfill \\ =gM\overrightarrow{X}·\hat{z},\hfill & \left(2.85\right)\hfill \end{array}$$

where the last sum is zero because the center of mass is the origin of ${\overrightarrow{\xi }}_{\alpha }$. So the potential energy of a body in a gravitational field with uniform acceleration is very simple: it is just M gh, where M is the total mass and $h=\overrightarrow{X}·\hat{z}$ is the height of the center of mass.

Here we consider an axisymmetric top (see figure 2.4). Such a top has an axis of symmetry of the mass distribution, so the center of mass is on the symmetry axis and the fixed point is also on the axis of symmetry.

In order to write the Lagrangian we need to choose a set of generalized coordinates. If we choose them well we can take advantage of the symmetries of the problem. If the Lagrangian does not depend on a particular coordinate, the conjugate momentum is conserved, and the complexity of the system is reduced.

The axisymmetric top has two apparent symmetries. The fact that the mass distribution is axisymmetric implies that neither the kinetic nor the potential energy is sensitive to the orientation of the top about that symmetry axis. Additionally, the kinetic and potential energy are insensitive to a rotation of the physical system about the vertical axis, because the gravitational field is uniform.

We can take advantage of these symmetries by choosing appropriate coordinates, and we already have a coordinate system that does the job—the Euler angles.¹⁴ We choose the reference orientation so that the symmetry axis is vertical. The first Euler angle, ψ, expresses a rotation about the symmetry axis. The next Euler angle, θ, is the tilt of the symmetry axis of the top from the vertical. The third Euler angle, φ, expresses a rotation of the top about the ẑ axis. The symmetries of the problem imply that the first and third Euler angles do not appear in the Lagrangian. As a consequence the momenta conjugate to these angles are conserved quantities. Let’s work out the details.

First, we develop the Lagrangian explicitly. The general form of the kinetic energy about a fixed point is given by equation 2.41. The top is constrained so that it pivots about a fixed point that is not at the center of mass. So the moments of inertia that enter the kinetic energy are the moments of inertia of the top with respect to the pivot point, not with respect to the center of mass. If we know the moments of inertia about the center of mass we can write the moments of inertia about the pivot in terms of them (see exercise 2.2 on Steiner’s theorem). So let’s assume the principal moments of inertia of the top about the pivot are A, B, and C, and A = B because of the symmetry.¹⁵ We can use the computer to help us figure out the Lagrangian for this special case:

We can rearrange this a bit to get

$$\begin{array}{ll}T\left(t;\theta ,\phi ,\psi ;\dot{\theta },\dot{\phi },\dot{\psi }\right)\hfill & \hfill \\ =\frac{1}{2}A\left({\dot{\theta }}^2+{\dot{\phi }}^2{\sin }^2\theta \right)+\frac{1}{2}C{\left(\dot{\psi }+\dot{\phi }\cos \theta \right)}^2.\hfill & \left(2.86\right)\hfill \end{array}$$

In terms of Euler angles, the potential energy is

$$\begin{array}{cc}V\left(t;\theta ,\phi ,\psi ;\dot{\theta },\dot{\phi },\dot{\psi }\right)=MgR\cos \theta ,& \left(2.87\right)\end{array}$$

where R is the distance of the center of mass from the pivot. The Lagrangian is L = T − V. We see that the Lagrangian is indeed independent of ψ and φ, as expected.

There is no particular reason to look at the Lagrange equations. We can assign that job to the computer when needed. However, we have already seen that it may be useful to examine the conserved quantities associated with the symmetries.

The energy is conserved, because the Lagrangian has no explicit time dependence. Also, the energy is the sum of the kinetic and potential energy E = T + V, because the kinetic energy is a homogeneous quadratic form in the generalized velocities. The energy is

$$\begin{array}{cc}E=\frac{1}{2}A\left({\dot{\theta }}^2+{\dot{\phi }}^2{\sin }^2\theta \right)+\frac{1}{2}C{\left(\dot{\psi }+\dot{\phi }\cos \theta \right)}^2+MgR\cos \theta .& \left(2.88\right)\end{array}$$

Two of the generalized coordinates do not appear in the Lagrangian, so there are two conserved momenta. The momentum conjugate to φ is

$$\begin{array}{cc}{p}_{\phi }=\left(A{\left(\sin \theta \right)}^2+C{\left(\cos \theta \right)}^2\right)\dot{\phi }+C\dot{\psi }\cos \theta .& \left(2.89\right)\end{array}$$

The momentum conjugate to ψ is

$$\begin{array}{cc}{p}_{\psi }=C\left(\dot{\psi }+\dot{\phi }\cos \theta \right).& \left(2.90\right)\end{array}$$

The state of the system at a moment is specified by the tuple $\left(t;\theta ,\phi ,\psi ;\dot{\theta },\dot{\phi },\dot{\psi }\right)$. Because the two coordinates φ and ψ do not appear in the Lagrangian, they do not appear in the Lagrange equations or the conserved momenta. So the evolution of the remaining four state variables, θ, $\dot{\theta }$, $\dot{\phi }$, and $\dot{\psi }$, depends only on those remaining state variables. This subsystem for the top has a four-dimensional state space. The variables that did not appear in the Lagrangian can be determined by integrating the derivatives of these variables, which are determined separately by solving the independent subsystem.

The evolution of the top is described by a four-dimensional subsystem and two auxiliary quadratures.¹⁶ This subdivision is a consequence of choosing generalized coordinates that incorporate the symmetries. However, the choice of generalized coordinates that incorporate the symmetries also gives conserved momenta. We can make use of these momenta to simplify the formulation of the problem further. Each conserved quantity can be used to locally eliminate one dimension of the subsystem. In this case the subsystem has four dimensions and there are three conserved quantities, so the system can be completely reduced to quadratures. For the top, this can be done analytically, but we think it is a waste of time to do so. Rather, we are interested in extracting interesting features of the motion. We concentrate on the energy and use the two conserved momenta to eliminate $\dot{\phi }$ and $\dot{\psi }$. After a bit of algebra we find:

$$\begin{array}{cc}E=\frac{1}{2}A{\dot{\theta }}^2+\frac{{\left(p_{\phi }-p_{\psi }\cos \theta \right)}^2}{2A{\left(\sin \theta \right)}^2}+\frac{p_{\psi }^2}{2C}+MgR\cos \theta .& \left(2.91\right)\end{array}$$

Along a path θ, where Dθ(t) is substituted for $\dot{\theta }$, this is an ordinary differential equation for θ. This differential equation involves various constants, some of which are set by the initial conditions of the other state variables. The solution of the differential equation for θ involves no more than ordinary integrals. So the top is essentially solved. We could continue this argument to obtain the qualitative behavior of θ: Using the energy (2.91), we can plot the trajectories in the plane of $\dot{\theta }$ versus θ and see that the motion of θ is simply periodic. However, we will defer this until chapter 3, when we have developed more tools for analysis.

Let’s get real. Let’s make a top out of an aluminum disk with a steel rod through the center to make the pivot. Measuring the top very carefully, we find that the moment of inertia of the top about the symmetry axis is about 1.32 × 10⁻⁴ kg m², and the moment of inertia about the pivot point is about 6.96 × 10⁻⁴ kg m². The combination gM R is about 0.112kg m² s⁻². We spin the top up with an initial angular velocity of $\dot{\psi }$ = 200 rad s⁻¹ (about 1910 rpm).

The top initially has $\dot{\theta }$ = φ = ψ = 0 and is initially tilted with θ = 0.4 rad. We then kick it so that $\dot{\phi }$ = −10 rad s⁻¹. Figures 2.5–2.8 display aspects of the evolution of the top for 2 seconds. The tilt of the top (measured by θ) varies in a periodic manner. The orientation about the vertical is measured by φ: we see that the top also precesses, and the rate of precession varies with θ. We also see that as the top bobs up and down the rate of rotation of the top oscillates—the top spins faster when it is more vertical. The plot of tilt versus precession angle shows that in this case the top executes a looping motion. If we do not kick it but just let it drop, then the loop disappears, leaving just a cusp. If we kick it in the other direction, then there is no cusp nor any looping motion.

Exercise 2.14: Kinetic energy of the top

The rotational kinetic energy of the top can be written in terms of the principal moments of inertia with respect to the pivot point and the angular velocity vector of rotation with respect to the pivot point. Show that this formulation of the kinetic energy yields the same value that one would obtain by computing the sum of the rotational kinetic energy about its center of mass and the kinetic energy of the motion of the center of mass.

Exercise 2.15: Nutation of the top

a. Carry out the algebra to obtain the energy (2.91) in terms of θ and $\dot{\theta }$.

b. Numerically integrate the Lagrange equations for the top to obtain figure 2.5, θ versus time.

c. Note that the energy is a differential equation for $\dot{\theta }$ in terms of θ, with conserved quantities p_φ, p_ψ, and E determined by initial conditions. Can we use this differential equation to obtain θ as a function of time? Explain.

Exercise 2.16: Precession of the top

Consider a top that is rotating so that θ is constant.

a. Using conservation of angular momentum, compute the rate of precession $\dot{\phi }$ as a function of the conserved angular momenta and the equilibrium value of θ.

b. For θ to be at an equilibrium the acceleration D²θ must be zero. Use the Lagrange equation for θ to find the rate of precession $\dot{\phi }$ at the equilibrium in terms of the equilibrium θ and $\dot{\psi }$.

c. Find an approximate expression for the precession rate in the limit that $\dot{\psi }$ is large.

d. The Newtonian rule is that the rate of change of the angular momentum is the torque. Assume the top is spinning so fast that the angular momentum is nearly the same as the angular momentum of rotation about the symmetry axis. By equating the rate of change of this vector angular momentum to the gravitational torque on the center of mass develop an approximate formula for the precession rate.

e. Numerically integrate the top to check your deductions.

2.11   Spin-Orbit Coupling

The rotation of planets and natural satellites is affected by the gravitational forces from other celestial bodies. As an extended application of the Lagrangian method for forced rigid bodies, we consider the rotation of celestial objects subject to gravitational forces.

We first develop the form of the potential energy for the gravitational interaction of an extended body with an external point mass. With this potential energy and the usual rigid-body kinetic energy we can form Lagrangians that model a number of systems. We will take an initial look at the rotation of the Moon and Mercury; later, after we have developed more tools, we will return to study these systems.

2.11.1 Development of the Potential Energy

The first task is to develop convenient expressions for the gravitational potential energy of the interaction of a rigid body with a distant point mass. A rigid body can be thought of as made of a large number of mass elements, subject to rigid coordinate constraints. We have seen that the kinetic energy of a rigid body is conveniently expressed in terms of the moments of inertia of the body and the angular velocity vector, which in turn can be represented in terms of a suitable set of generalized coordinates. The potential energy can be developed in a similar manner. We first represent the potential energy in terms of moments of the mass distribution and later introduce generalized coordinates as particular parameters of the potential energy.

The gravitational potential energy of a point mass and a rigid body (see figure 2.9) is the sum of the potential energy of the point mass with each mass element of the body:

$$\begin{array}{cc}-{\displaystyle \sum _{\alpha }\frac{GM\prime m_{\alpha }}{r_{\alpha }}},& \left(2.92\right)\end{array}$$

where M′ is the mass of the external point mass, r_α is the distance between the point mass and the constituent mass element with index α, m_α is the mass of this constituent element, and G is the gravitational constant. Let R be the distance of the center of mass of the rigid body from the point mass. The distance from the center of mass to the constituent with index α is ξ_α. The distance r_α is then given by the law of cosines as $r_{\alpha }^2=R^2+{\xi }_{\alpha }^2-2{\xi }_{\alpha }R\cos {\theta }_{\alpha }$, where θ_α is the angle between the lines from the center of mass to the constituent and to the point mass.

Because this is a three-dimensional body the distance ξ_α and angle θ_α do not completely specify the position of the constituent mass element; to do that one must also specify the angle of rotation about the line between the center of mass and the external point mass. But the potential energy does not depend on this angle.

The potential energy is then

$$\begin{array}{cc}-GM\prime {\displaystyle \sum _{\alpha }\frac{m_{\alpha }}{{\left(R^2+{\xi }_{\alpha }^2-2{\xi }_{\alpha }R\cos {\theta }_{\alpha }\right)}^{1/2}}.}& \left(2.93\right)\end{array}$$

This is complete, but we need to find a representation that does not mention each constituent.

Typically, the size of celestial bodies is small compared to the distance between them. We can make use of this to find a more compact representation of the potential energy. If we expand the potential energy in the small ratio ξ_α/R we find

$$\begin{array}{cc}-GM\prime {\displaystyle \sum _{\alpha }{m}_{\alpha }\frac{1}{R}{\displaystyle \sum _l\frac{{\xi }_{\alpha }^l}{R^l}{P}_l\left(\cos {\theta }_{\alpha }\right),}}& \left(2.94\right)\end{array}$$

where P_l is the lth Legendre polynomial.¹⁷ Interchanging the order of the summations yields:

$$\begin{array}{cc}-\frac{GM\prime }{R}{\displaystyle \sum _l{\displaystyle \sum _{\alpha }{m}_{\alpha }}\frac{{\xi }_{\alpha }^l}{R^l}{P}_l\left(\cos {\theta }_{\alpha }\right).}& \left(2.95\right)\end{array}$$

Successive terms in this expansion of the potential energy typically decrease very rapidly because celestial bodies are small compared to the distance between them. We can compute an upper bound to the size of these terms by replacing each factor in the sum over α by an upper bound. The Legendre polynomials all have magnitudes less than one for arguments in the range −1 to 1. The distances ξ_α are all less than some maximum extent of the body ξ_max. The sum over m_α times these upper bounds is just the total mass M times the upper bounds. Thus

$$\begin{array}{cc}\left|{\displaystyle \sum _{\alpha }{m}_{\alpha }\frac{{\xi }_{\alpha }^l}{R^l}{P}_l\left(\cos {\theta }_{\alpha }\right)}\right|\le M\frac{{\xi }_{\max }^l}{R^l}.& \left(2.96\right)\end{array}$$

We see that the upper bound on successive terms decreases by a factor ξ_max/R. Successive terms may be smaller still. For large bodies the gravitational force is strong enough to overcome the internal material strength of the body, so the body, over time, becomes nearly spherical. Successive terms in the expansion of the potential are measures of the deviation of the mass distribution from a spherical mass distribution. Thus for large bodies the higher-order terms are small because the bodies are nearly spherical.

Consider the first few terms in l. For l = 0 the sum over α just gives the total mass M of the rigid body. For l = 1 the sum over α is zero, as a consequence of choosing the origin of the ${\overrightarrow{\xi }}_{\alpha }$ to be the center of mass. For l = 2 we have to do a little more work. The sum involves second moments of the mass distribution, and can be written in terms of moments of inertia of the rigid body:

$$\begin{array}{lll}{\displaystyle \sum _{\alpha }{m}_{\alpha }{\xi }_{\alpha }^2{P}_2\left(\cos {\theta }_{\alpha }\right)}\hfill & ={\displaystyle \sum _{\alpha }{m}_{\alpha }{\xi }_{\alpha }^2\left(\frac{3}{2}{\left(\cos {\theta }_{\alpha }\right)}^2-\frac{1}{2}\right)}\hfill & \hfill \\ \hfill & ={\displaystyle \sum _{\alpha }{m}_{\alpha }{\xi }_{\alpha }^2\left(1-\frac{3}{2}{\left(\sin {\theta }_{\alpha }\right)}^2\right)}\hfill & \hfill \\ \hfill & =\frac{1}{2}\left(A+B+C-3I\right),\hfill & \left(2.97\right)\hfill \end{array}$$

where A, B, and C are the principal moments of inertia, and I is the moment of inertia of the rigid body about the line between the center of mass of the body and the external point mass. The moment I depends on the orientation of the rigid body relative to the line between the bodies. The contributions to the potential energy up to l = 2 are then¹⁸

$$\begin{array}{cc}-\frac{GMM\prime }{R}-\frac{GM\prime }{2R^3}\left(A+B+C-3I\right).& \left(2.98\right)\end{array}$$

Let c_a = cos θ_a, c_b = cos θ_b, and c_c = cos θ_c be the direction cosines of the angles θ_a, θ_b and θ_c between the principal axes â, $\hat{b}$, and ĉ and the line between the center of mass and the point mass. (See figure 2.10.) A little algebra shows that $I=c_a^{2}A+c_b^{2}B+c_c^{2}C$. The potential energy is then

$$\begin{array}{cc}-\frac{GMM\prime }{R}-\frac{GM\prime }{2R^3}\left[\left(1-3c_a^2\right)A+\left(1-3c_b^2\right)B+\left(1-3c_c^2\right)C\right].& \left(2.99\right)\end{array}$$

This is a good first approximation to the potential energy of interaction for most situations in the solar system; if we intended to land on the Moon we probably would want to take into account higher-order terms in the expansion.

Exercise 2.17:

a. Fill in the details that show that the sum over constituents in equation (2.97) can be expressed as written in terms of moments of inertia. In particular, show that

$\begin{array}{l}{\displaystyle \sum _{\alpha }{m}_{\alpha }{\xi }_{\alpha }\cos {\theta }_{\alpha }=0,}\\ 2{\displaystyle \sum _{\alpha }{m}_{\alpha }{\xi }_{\alpha }^2=A+B+C,}\end{array}$

and that

${\displaystyle \sum _{\alpha }{m}_{\alpha }{\xi }_{\alpha }^2{\left(\sin {\theta }_{\alpha }\right)}^2=I.}$

b. Show that if the principal moments of inertia of a rigid body are A, B, and C, then the moment of inertia about an axis that goes through the center of mass of the body with angles θ_a, θ_b, and θ_c to the principal axes is

$I={\left(\cos {\theta }_a\right)}^{2}A+{\left(\cos {\theta }_b\right)}^{2}B+{\left(\cos {\theta }_c\right)}^{2}C.$

2.11.2 Rotation of the Moon and Hyperion

The approximation to the potential energy that we have derived can be used for a number of different problems. It can be used to investigate the effect of oblateness on the motion of an artificial satellite about the Earth, or to incorporate the effect of planetary oblateness on the evolution of the orbits of natural satellites, such as the Moon or the Galilean satellites of Jupiter. However, as the principal application here, we will use it to investigate the rotational dynamics of natural satellites and planets.

The potential energy depends on the position of the point mass relative to the rigid body and on the orientation of the rigid body. Thus the changing orientation is coupled to the orbital evolution; each affects the other. However, in many situations the effect of the orientation of the body on the evolution of the orbit may be ignored. One way to see this is to look at the relative magnitudes of the two terms in the potential energy (2.99). We already know that the second term is guaranteed to be smaller than the first by a factor of (ξ_max/R)², but often it is much smaller still because the body involved is nearly spherical. For example, the radius of the Moon is about a third the radius of the Earth and the distance to the Moon is about 60 Earth-radii. So the second term is smaller than the first by a factor of order 10⁻⁴ due to the size factors. In addition, the Moon is roughly spherical and for any orientation the combination A + B + C − 3I is of order 10⁻⁴C. Now C is itself of order $\frac{2}{5}M{R}^2$, because the density of the Moon does not vary strongly with radius. So for the Moon the second term is of order 10⁻⁸ relative to the first. Even radical changes in the orientation of the Moon would have little dynamical effect on its orbit.

We can learn some important qualitative aspects of the orientation dynamics by studying a simplified model problem. First, we assume that the body is rotating about its largest moment of inertia. This is a natural assumption. Remember that for a free rigid body the loss of energy while conserving angular momentum leads to rotation about the largest moment of inertia. This is observed for most bodies in the solar system. Next, we assume that the spin axis is perpendicular to the orbital motion. This is a good approximation for the rotation of natural satellites, and is a natural consequence of tidal friction—dissipative solid-body tides raised on the satellite by the gravitational interaction with the planet. Finally, for simplicity we take the rigid body to be moving on a fixed elliptic orbit. This may approximate the motion of some physical systems, provided the time scale of the evolution of the orbit is large compared to any time scale associated with the rotational dynamics we are investigating. So we have a nice toy problem, one that has been used to investigate the rotational dynamics of Mercury, the Moon, and other natural satellites. It makes specific predictions concerning the rotation of Phobos, a satellite of Mars, that can be compared with observations. It provides a basic explanation of the fact that Mercury rotates precisely three times for every two orbits it completes, and is the starting point for understanding the chaotic tumbling of Saturn’s satellite Hyperion.

We are assuming that the orbit does not change or precess. The orbit is an ellipse with the point mass at a focus of the ellipse. The angle f (see figure 2.11) measures the position of the rigid body in its orbit relative to the point in the orbit at which the two bodies are closest.¹⁹ We assume the orbit is a fixed ellipse, so the angle f and the distance R are periodic functions of time, with period equal to the orbit period. With the spin axis constrained to be perpendicular to the orbit plane, the orientation of the rigid body is specified by a single degree of freedom: the orientation of the body about the spin axis. We specify this orientation by the generalized coordinate θ that measures the angle to the â principal axis from the same line from which we measure f, the line through the point of closest approach.

Having specified the coordinate system, we can work out the details of the kinetic and potential energies, and thus find the Lagrangian. The kinetic energy is

$$\begin{array}{cc}T\left(t,\theta ,\dot{\theta }\right)=\frac{1}{2}C{\dot{\theta }}^2,& \left(2.100\right)\end{array}$$

where C is the moment of inertia about the spin axis and the angular velocity of the body about the ĉ axis is $\dot{\theta }$. There is no component of angular velocity on the other principal axes.

To get an explicit expression for the potential energy, write the direction cosines in terms of θ and f: cos θ_a = − cos(θ − f), cos θ_b = sin(θ − f), and cos θ_c = 0 because the ĉ axis is perpendicular to the orbit plane. The potential energy is then

$$\begin{array}{l}-\frac{GMM\prime }{R}\hfill \\ -\frac{1}{2}\frac{GM\prime }{R^3}\left[\left(1-3{\cos }^2\left(\theta -f\right)\right)A+\left(1-3{\sin }^2\left(\theta -f\right)\right)B+C\right].\hfill \end{array}$$

Since we are assuming that the orbit is given, we need keep only terms that depend on θ. Expanding the squares of the cosine and the sine in terms of the double angles and dropping all the terms that do not depend on θ, we find the potential energy for the orientation²⁰

$$\begin{array}{cc}V\left(t,\theta ,\dot{\theta }\right)=-\frac{3}{4}\frac{GM\prime }{R^3\left(t\right)}\left(B-A\right)\cos 2\left(\theta -f\left(t\right)\right).& \left(2.101\right)\end{array}$$

A Lagrangian for the model spin-orbit coupling problem is then L = T − V :

$$\begin{array}{cc}L\left(t,\theta ,\dot{\theta }\right)=\frac{1}{2}C{\dot{\theta }}^2+\frac{3}{4}\frac{GM\prime }{R^3\left(t\right)}\left(B-A\right)\cos 2\left(\theta -f\left(t\right)\right).& \left(2.102\right)\end{array}$$

We introduce the dimensionless “out-of-roundness” parameter

$$\begin{array}{cc}\mathit{ϵ}=\sqrt{\frac{3\left(B-A\right)}{C},}& \left(2.103\right)\end{array}$$

and use the fact that the orbital frequency n and the semimajor axis a satisfy Kepler’s third law, n²a³ = G(M + M′), which is approximately n²a³ = GM′ for a small body in orbit around a much more massive one (M ≪ M′). In terms of ϵ and n the spin-orbit Lagrangian is

$$\begin{array}{cc}L\left(t,\theta ,\dot{\theta }\right)=\frac{1}{2}C{\dot{\theta }}^2+\frac{n^2{\mathit{ϵ}}^{2}C}{4}\frac{a^3}{R^3\left(t\right)}\cos 2\left(\theta -f\left(t\right)\right).& \left(2.104\right)\end{array}$$

This is a problem with one degree of freedom with terms that vary periodically with time.

The Lagrange equations are derived in the usual manner:

$$\begin{array}{cc}C{D}^2\theta \left(t\right)=-\frac{n^2{\mathit{ϵ}}^{2}C}{2}\frac{a^3}{R^3\left(t\right)}\sin 2\left(\theta \left(t\right)-f\left(t\right)\right).& \left(2.105\right)\end{array}$$

The equation of motion is very similar to that of the periodically driven pendulum. The main difference here is that not only is the strength of the acceleration changing periodically, but in the spin-orbit problem the center of attraction is also varying periodically.

We can give a physical interpretation of this equation of motion. It states that the rate of change of the angular momentum is equal to the applied torque. The torque on the body arises because the body is out of round and the gravitational force varies as the inverse square of the distance. Thus the force per unit mass on the near side of the body is a little more than the acceleration of the body as a whole, and the force per unit mass on the far side of the body is a little less than the acceleration of the body as a whole. Thus, relative to the acceleration of the body as a whole, the far side is forced outward while the inner part of the body is forced inward. The net effect is a torque on the body that tries to align the long axis of the body with the line to the external point mass. If θ is a bit larger than f then there is a negative torque, and if θ is a bit smaller than f then there is a positive torque, both of which would align the long axis with the point mass if given a fair chance. The torque arises because of the difference of the inverse R² force across the body, so the torque is proportional to R⁻³. There is a torque only if the body is out of round, for otherwise there is no handle to pull on. This is reflected in the factor B − A in the expression for the torque. The potential depends on the mass distribution as described by the moments of inertia, and thus the body has the same dynamics if it is rotated by 180°. The factor of 2 in the argument of sine reflects this symmetry. This torque is called the “gravity gradient torque.”

To compute the evolution requires a lot of detailed preparation similar to what has been done for other problems. There are many interesting phenomena to explore. We can take parameters appropriate for the Moon and find that Mr. Moon does not steadily point the same face to the Earth, but instead constantly shakes his head in dismay at what goes on here. If we nudge the Moon a bit, say by hitting it with an asteroid, we find that the long axis oscillates back and forth with respect to the direction that points to the Earth. For the Moon, the orbital eccentricity is currently about 0.05, and the out-of-roundness parameter is about ϵ = 0.026. Figure 2.12 shows the angle θ − f as a function of time for two different values of the “lunar” eccentricity. The plot spans 50 lunar orbits, or a little under four years. This Moon has been kicked by a large asteroid and has initial rotational angular velocity $\dot{\theta }$ equal to 1.01 times the orbit frequency. The initial orientation is θ = 0. The smooth trace shows the evolution if the orbital eccentricity is set to zero. We see an oscillation with a period of about 40 lunar orbits or about three years. The more wiggly trace shows the evolution of θ − f with an orbital eccentricity of 0.05, near the current lunar eccentricity. The lunar eccentricity superimposes an apparent shaking of the face of the Moon back and forth with the period of the lunar orbit. Though the Moon does slightly change its rate of rotation during the course of its orbit, most of this shaking is due to the nonuniform motion of the Moon in its elliptical orbit. This oscillation, called the “optical libration of the Moon,” allows us to see a bit more than half of the Moon’s surface. The longer-period oscillation induced by the kick is called the “free libration of the Moon.” It is “free” because we are free to excite it by choosing appropriate initial conditions. The mismatch of the orientation of the Moon caused by the optical libration actually produces a periodic torque on the Moon, which slightly speeds it up and slows it down during every orbit. The resulting oscillation is called the “forced libration of the Moon,” but it is too small to see in this plot.

The oscillation period of the free libration is easily calculated. We see that the eccentricity of the orbit does not substantially affect the period, so we consider the special case of zero eccentricity. In this case R = a, a constant, and f(t) = nt, where n is the orbital frequency.²¹ The equation of motion becomes

$$\begin{array}{cc}{D}^2\theta \left(t\right)=-\frac{n^2{\mathit{ϵ}}^2}{2}\sin 2\left(\theta \left(t\right)-nt\right).& \left(2.106\right)\end{array}$$

Let φ(t) = θ(t) − nt, and consequently Dφ(t) = Dθ(t) − n, and D²φ = D²θ. Substituting these, the equation governing the evolution of φ is

$$\begin{array}{cc}{D}^2\phi =-\frac{n^2{\mathit{ϵ}}^2}{2}\sin 2\phi .& \left(2.107\right)\end{array}$$

For small deviations from synchronous rotation (small φ) this is

$$\begin{array}{cc}{D}^2\phi =-n^2{\mathit{ϵ}}^2\phi ,& \left(2.108\right)\end{array}$$

so we see that the small-amplitude oscillation frequency of φ is nϵ. For the Moon, ϵ is about 0.026, so the period is about 1/0.026 orbit periods or about 40 lunar orbit periods, which is what we observed.

It is perhaps more fun to see what happens if the out-of-roundness parameter is large. After our experience with the driven pendulum it is no surprise that we find abundant chaos in the spin-orbit problem when the system is strongly driven by having large ϵ and significant orbital eccentricity e. There is indeed one body in the solar system that exhibits chaotic rotation—Hyperion, a small satellite of Saturn. Though our toy model is not adequate for a complete account of Hyperion, we can show that it exhibits chaotic behavior for parameters appropriate for Hyperion. We take ϵ = 0.89 and e = 0.1. Figure 2.13 shows θ − f for 50 orbits, starting with θ = 0 and $\dot{\theta }$ = 1.05. We see that sometimes one face of the body oscillates facing the planet, sometimes the other face oscillates facing the planet, and sometimes the body rotates relative to the planet in either direction.

If we relax our restriction that the spin axis be fixed perpendicular to the orbit, then we find that the Moon maintains this orientation of the spin axis even if nudged a bit, but for Hyperion the spin axis almost immediately falls away from this configuration. The state in which Hyperion on average points one face to Saturn is dynamically unstable to chaotic tumbling. Observations of Hyperion are consistent with the deduction that it is chaotically tumbling.

Exercise 2.18: Precession of the equinox

The Earth spins very nearly about the largest moment of inertia, and the spin axis is tilted by about 23° to the orbit normal. There is a gravity-gradient torque on the Earth from the Sun that causes the spin axis of the Earth to precess. Investigate this precession in the approximation that the orbit of the Earth is circular and the Earth is axisymmetric. Determine the rate of precession in terms of the moments of inertia of the Earth.

2.11.3 Spin-Orbit Resonances

Consider the motion of the Moon in synchronous rotation. We have seen that if we give the Moon a kick so that it is not exactly pointing one face to the Earth, then the face will oscillate back and forth relative to the direction to the Earth. If we give it a really big kick, then instead of oscillating it will spin relative to the direction to the Earth. How do we understand this?

Let’s look again at the equations of motion for the rotation of the Moon when the orbit is circular (equation 2.106):

$$\begin{array}{cc}C{D}^2\theta (t)=-\frac{C}{2}{n}^2{\mathit{ϵ}}^2\sin 2(\theta (t)-n)t).& \left(2.109\right)\end{array}$$

Changing variables to φ(t) = θ(t) − nt this equation becomes

$$\begin{array}{cc}C{D}^2\phi \left(t\right)=-\frac{C}{2}{n}^2{\mathit{ϵ}}^2\sin 2\phi \left(t\right).& \left(2.110\right)\end{array}$$

This equation can be solved; it has an “energy-like” conserved quantity

$$\begin{array}{cc}E\left(\phi ,\dot{\phi }\right)=\frac{C}{2}{\dot{\phi }}^2-\frac{C}{4}{n}^2{\mathit{ϵ}}^2\cos \left(2\phi \right).& \left(2.111\right)\end{array}$$

The solutions are just contours of this conserved quantity (see figure 2.14). There are two centers of oscillation corresponding to the two different faces of the Moon that could point towards Earth. There are also trajectories that rotate relative to the Earth. And there are separating trajectories that divided the oscillating trajectories from the circulating trajectories. Where these separating trajectories appear to cross, the system is at an unstable equilibrium. The separating trajectories are asymptotic to the unstable equilibria (a system on that trajectory takes an infinite time to get to the equilibrium point). These asymptotic trajectories are analogous to the trajectories of the simple pendulum that are asymptotic to the vertical.

The extent of the oscillation region can be evaluated with the help of the conserved quantity E. Let’s evaluate it on the separating trajectory:

$$\begin{array}{lll}E\left(\pi /2,0\right)\hfill & =\frac{C}{4}{n}^2{\mathit{ϵ}}^2,\hfill & \hfill \\ E\left(0,\Delta \dot{\phi }\right)\hfill & =\frac{C}{2}\Delta {\dot{\phi }}^2-\frac{C}{4}{n}^2{\mathit{ϵ}}^2.\hfill & \left(2.112\right)\hfill \end{array}$$

Equating these and solving, we find the maximum extent of the oscillating region:

$$\begin{array}{cc}\Delta \dot{\phi }=n\mathit{ϵ}.& \left(2.113\right)\end{array}$$

So we see that the out-of-roundness parameter ϵ not only gives the frequency of small-amplitude oscillations, but also gives the extent of the oscillation region.

Mercury rotates exactly three times for every two times it goes around the Sun, as discovered by Pettengill and Dyce in 1965, using Arecibo radar. We can understand this spin-orbit resonance using our simple spin-orbit model problem.

Let’s first use qualitative reasoning to understand how the resonance comes about. The spin-orbit equation of motion, equation (2.105), equates the rate of change of the angular momentum to the gravity gradient torque. The torque is proportional to the inverse cube of the distance, so it is largest when the distance is smallest, at pericenter. For the purpose of qualitative reasoning, consider the torque only at pericenter. Suppose Mercury is rotating exactly three times for every two orbits. Then if the long axis of Mercury is pointed at the Sun at one pericenter, it will point the other end of this long axis the next time it passes pericenter (it will have rotated one and a half times). Now, suppose Mercury is rotating a little faster. Then if the long axis is aligned at one pericenter passage, on the next pericenter passage it will have rotated a little too much and the long axis will no longer point to the Sun. In this case θ − f will be positive and there will be a negative torque, slowing down the rotation a bit. Over many orbits the rotation of Mercury is reduced. A similar argument shows that if Mercury is rotating slower than three times for every two orbits, then the torques at succesive pericenter passages will tend to increase the rotation rate. An oscillation ensues.

We can also understand this spin-orbit resonance analytically. The right-hand side of the equation of motion (equation 2.105) has factors that vary periodically with the orbit period:

$$\begin{array}{cc}C{D}^2\theta \left(t\right)=-\frac{C}{2}{n}^2{\mathit{ϵ}}^2{\left(\frac{a}{R\left(t\right)}\right)}^3\sin 2\left(\theta \left(t\right)-f\left(t\right)\right),& \left(2.114\right)\end{array}$$

where both f(t) and R(t) are periodic with period 2π/n. We can expand this as a Fourier series:

$$\begin{array}{cc}C{D}^2\theta \left(t\right)=-\frac{C}{2}{n}^2{\mathit{ϵ}}^2{\displaystyle \sum _{m=-\infty }^{\infty }{A}_m\left(e\right)\sin \left(2\theta \left(t\right)-mnt\right),}& \left(2.115\right)\end{array}$$

where the coefficients A_m(e) are functions of the orbital eccentricity e. The coefficients are proportional to e^|m−2| and so for small eccentricity we need to consider only a few of them.²² We have

$$\begin{array}{lll}{A}_1(e)\hfill & =-\frac{e}{2}+o(e^3)\hfill & \hfill \\ A_2(e)\hfill & =1-\frac{5e^2}{2}+o(e^4)\hfill & \hfill \\ A_3(e)\hfill & =\frac{7e}{2}+o(e^3).\hfill & (2.116)\hfill \end{array}$$

All other terms are of higher order in e. With just the terms of order e or less, the equation of motion becomes

$$\begin{array}{lll}C{D}^2\theta \left(t\right)\hfill & =-\frac{C}{2}{n}^2{\mathit{ϵ}}^2[\sin \left(2\theta \left(t\right)-2nt\right)\hfill & \hfill \\ \hfill & +\frac{7e}{2}\sin \left(2\theta \left(t\right)-3nt\right)\hfill & \hfill \\ \hfill & -\frac{e}{2}\sin \left(2\theta \left(t\right)-nt\right)\hfill & \hfill \\ \hfill & +\cdots ].\hfill & \left(2.117\right)\hfill \end{array}$$

Suppose we are close to the 3:2 Mercury resonance. Then $\dot{\theta }$ is close to (3/2)n. So the combination θ(t) − (3/2)nt is slowly varying compared to the other two arguments: θ(t)−nt and θ(t)− (1/2)nt. The rapidly varying torques due to these other terms tend to average out, leaving a slowly varying torque that controls the motion.²³ The averaged equation of motion for motion near the 3:2 resonance is then

$$\begin{array}{cc}C{D}^2\theta \left(t\right)=-\frac{C}{2}{n}^2{\mathit{ϵ}}^2\left(\frac{7e}{2}\right)\sin (2\theta (t)-3nt).& \left(2.118\right)\end{array}$$

We can solve this by changing variables to

$$\begin{array}{cc}\phi \left(t\right)=\theta \left(t\right)-\left(3/2\right)nt.& \left(2.119\right)\end{array}$$

The equation of motion becomes

$$\begin{array}{cc}C{D}^2\phi \left(t\right)=\frac{C}{2}{n}^2{\mathit{ϵ}}^2\left(\frac{7e}{2}\right)\sin \left(2\phi \left(t\right)\right).& \left(2.120\right)\end{array}$$

This has the “energy-like” conserved quantity

$$\begin{array}{cc}E\left(\phi ,\dot{\phi }\right)=\frac{C}{2}{\dot{\phi }}^2-\frac{C}{4}{n}^2{\mathit{ϵ}}^2\left(\frac{7e}{2}\right)\cos \left(2\phi \right),& \left(2.121\right)\end{array}$$

which is very similar to the conserved quantity we found for the zero-eccentricity synchronous rotation case considered earlier; see equation (2.111). Indeed the trajectories are contours of the conserved quantity and look just like those in figure 2.14. Using analogous reasoning we can determine the extent of the oscillation region and find

$$\begin{array}{cc}\Delta \dot{\phi }=n\mathit{ϵ}\sqrt{\frac{7e}{2}.}& \left(2.122\right)\end{array}$$

This gives the approximate range of rotation rate over which Mercury can oscillate stably in the 3:2 resonance.