Appendix I on rigid body

Appendix I: Rigid Body Dynamics 1 1.1 The Appearance of the Inertia Tensor 1 I.2 Rigid Body Equations of Motion : Part 1 5 I.3 Diagonalization of the Inertia Tensor in Frame S' 6 I.4 Rigid Body Equations of Motion : Part 2 7 I.5 Zero-torque motion of a Rigid Body : Ellipsoids and Poinsot 8 I.6 Zero-torque motion of an Axisymmetric Rigid Body (Rigid Rotor) 15 I.7 Rigid Body with External Torque: Spinning Symmetric Top 25 I.8 Gravitational Torque on an oblate Earth: Precession of the Equinoxes 31 I.9 Derivation of the Oblate Earth torque formula 34 This has been installed, do not edit here. Appendix I: Rigid Body Dynamics We shall develop here the equations of motion for a rigid body's angular momentum vector ω(t). This is done using ω components (ω)'i in Frame S', the body frame. It is shown that one might as well choose the Frame S' axes to line up with some set of principal axes of the rigid body since this diagonalizes the inertia tensor, simplifying equations. We then examine the general torque-free solution for ω using the construction of Poinsot, and then we analytically solve for ω for the special axisymmetric case. The results for ω are obtained in both Frame S' and inertial Frame S and appropriate "cone pictures" are displayed. It is noted that the Earth exhibits a tiny torque-free precession called the Chandler wobble. We then consider the torque-present problem of a spinning top. For the Earth, weak torques are applied by the Sun and the Moon causing a slow precession known as the precession of the equinoxes. In the final section we derive an expression for the torque of the Sun or Moon on the slightly oblate Earth. 1.1 The Appearance of the Inertia Tensor Assume that Frame S' is embedded into a rigid body, so that Frame S' and that rigid body are rotating at angular velocity ω relative to an inertial Frame S. This is the situation of our "non-swap notation". Goldstein uses the appropriate term "body frame" for such a Frame S'. All sensible authors use "swap notation" for this analysis to avoid an avalanche of prime symbols (or other labels), but we shall be obstinate in using the non-swap notation because it forces us to think carefully about many details. We have to decide whether to put a prime or not put a prime on each entity. In Frame S' the rigid body is at rest, so the Frame S' velocity v'α of any particle α of the rigid body vanishes, v'α = 0. For particle α we also have p'α = mαv'α = 0 and L'α = r'α x p'α = 0. In particular, the total Frames S' angular momentum of the rigid body is L' = ΣαL'α = 0. These facts are quite obvious but we state them anyway. Note that α here is a label and not a component index. We shall use α in this manner as a particle label, and the usual i,j,k indices as component indices. Obviously, for a continuous rigid body Σα is really an integration over the particles of the body. In Frame S, however, the rigid body has some non-vanishing angular momentum L = Σα rα x mvα . In what follows, we shall be very careful with the use of "primes" to clearly show what objects are the natural objects in what frames, and what frame components are being evaluated. As noted in Section 1.3, in order to mark our "natural" objects with primes or no primes depending on whether they are associated with Frame S' or Frame S, we must use the Passive View of rotation transformations. For vectors this means (V)' = RV and for rank-2 tensors (T)' = RTR-1 (with Frame S' components (V)'i and (T)'ij ). Angular Momentum, the Inertia Tensor, and Kinetic Energy If we assume that Frame S and Frame S' have a common origin, we are in Special Case #4 discussed in Section 10 and summarized below Fig (12.1.7). In Special Case #4 the vector b of Fig 1 is 0, so r' = r. In this case, for a rigid body particle α one has, from (12.1.8b), vα = v'α + ω x r'α = 0 + ω x r'α = ω x rα . (12.1.8b) (I.1.1) Then the total Frame S angular momentum of the rigid body is L = Σα rα x pα = Σα rα x mαvα = Σα mα rα x (ω x rα) = Σα mα [ rα2ω - (rα ω) rα ] . //A x (B x C) = (AC)B - (AB)C (I.1.2) Taking Frame S components, Li = Σα mα [ rα2ωi - (rα ω) (rα)i ] = Σα mα [ rα2Σjωjδij - Σj(rα)jωj (rα)i ] = Σα mα Σj[ rα2δij - (rα)i(rα)j ] ωj = Σj { Σαmα [ rα2δij - (rα)i(rα)j ] } ωj = Σj Iijωj where Iij ≡ Σαmα [ rα2δij - (rα)i(rα)j ] = the inertia tensor (I.1.3) or Iij ≡ ∫dx1dx2dx3 ρ(x)[ r2δij - xixj] . // continuum notation We rewrite (I.1.3) as a vector equation, L = Iω . (I.1.4) Exercise: Write an expression for the Dirac inertia operator I (script I, not T) : <ei| I | ej> = Iij = Σαmα [ rα2δij - (rα)i(rα)j ] = Σαmα [ rα2<ei|ej> - <ei|rα><rα|ej> = <ei| { Σαmα [ rα2 1 - |rα><rα| } |ej> I = Σαmα ( rα21 - |rα><rα| ) . (I.1.5) In equation (I.1.4) L, I and ω are all associated with a Frame S observer. L is the "natural" angular momentum of the rigid body as seen in Frame S, and ω is the angular velocity of Frame S' (and its rigid body) as seen from Frame S. Finally Iij are the inertia tensor components viewed from Frame S. Although L = Iω is a Frame S equation, as with any vector equation one can take components of the equation in either Frame S or Frame S' : Li = Iij (ω)j Iij = Σαmα [ rα2δij - (rα)i(rα)j ] // rα2 = (rα)i(rα)i (I.1.6) (L)'i = (I)'ij (ω)'j (I)'ij = Σαmα [ r'α2δij - (rα)'i(rα)'j ] . // r'α2 = (rα)'i(rα)'i Notice in (I)'ij that (rα)'i = <e'i| rα > is constant in time. The components of a vector static in Frame S' do not change in time, although Frame S' moves relative to Frame S. Thus ∂t(rα)'i = 0. The magnitude of the vector rα is obviously constant in time, but we verify: ∂tr'α2 = ∂t[(rα)'i(rα)'i] = (rα)'∂t[(rα)'i] = 0. Therefore ∂t(I)'ij = 0 for the entire tensor (I'). In contrast, ∂tIij ≠ 0 . Exercise: Verify that (I)' = RIR-1, thus showing that I is a rank-2 tensor (in the Passive View) : (I)'ij = RiaIabR-1bj = Ria[ Σαmα ( rα2δab - (rα)a(rα)b )]R-1bj = Σαmα (rα2 RiaδabR-1bj - Ria(rα)aR-1bj(rα)b ) = Σαmα (rα2RiaR-1aj - [Ria(rα)a][Rjb(rα)b] ) = Σαmα (rα2 δij - (rα)'i(rα)'j) . QED In what follows, we ignore kinetic energy associated with motion of the rigid body center of mass. We think of this center of mass being at rest in both Frame S and Frame S' The total rotational Frame S kinetic energy T of the rigid body is given by, using (I.1.1), T = Σα(1/2)mαvα2 = Σα(1/2)mα vα (ω x rα) = Σα(1/2)mα ω (rα x vα) // cyclic rule = Σα (1/2) ω (rα x mαvα) = (1/2) ω Σα rα x pα = (1/2) ω L = (1/2) ω (Iω) . // (I.1.4) (I.1.7) We pause to note the various notations one can use for this dot product : ω (Iω) = ωTIω (matrix notation) = ωiIijωj = ωiωjIij (all Frame S components) = <ω | I | ω> (Dirac notation, where I is the abstract inertia operator = <ω |ei><ei| I |ej><ej| ω> = ωi Iij ωj // completeness (1.1.20) = ω I ω (Goldstein p 149 5-15) = ωω = ωIω . (I.1.8) Comments: We feel it is useful to keep stressing these notational issues: The Frame S' kinetic energy T' = (1/2) ω' I' ω' = 0 because ω' = 0, meaning (ω')i= 0. This is so because the rigid body is at rest in Frame S' -- nothing is moving in Frame S' ! The Frame S kinetic energy T = (1/2)ωIω can be written out in either Frame S or Frame S' components. T = (1/2) ωIω = (1/2) ωiωjIij = (1/2) (ω)'i(ω)'j(I)'ij . (I.1.9) Proof: In Dirac notation, T = (1/2) <ω | I | ω> = (1/2) <ω |ei><ei| I |ej><ej| ω> = (1/2) ωiIijωj T = (1/2) <ω | I | ω> = (1/2) <ω |e'i><e'i| I |e'j><e'j| ω> = (1/2)(ω)'i(I)'ij (ω)'j . Using the Goldstein double dot notation, kinetic energy can be rewritten as T = (1/2) ( I ) ω2 = (1/2) I ω2 where I ≡ I (I.1.10a) in analogy with T = (1/2)mv2 for linear motion. The scalar object I is called the moment of inertia of the rigid body about the axis, I ≡ I = ijIij = ij Σαmα [ rα2δij - (rα)i(rα)j ] = Σαmα [ rα2ijδij - (rα)i i(rα)jj ] = Σαmα [ rα2 - (rα )2] = Σαmα [ rα2 - (rα)||2] = Σαmα (rα)2, (I.1.10b) where (rα) is the perpendicular distance of particle α from the ω rotation axis, (I.1.10c) I.2 Rigid Body Equations of Motion : Part 1 From (11.3.2) Newton's Second Law (angular motion) applied in inertial Frame S states, N = // = ∂SL = (dL/dt)S (I.2.1) where we set the torque and angular momentum reference point to the common Frames' origin. This law is valid because Frame S is an inertial frame so there are no fictitious torques. N is the Frame S torque, and is the natural time derivative of natural L in Frame S (see (1.8.4) and (1.9.3)). Using the G Rule (2.1) one finds, in the notation of (1.8.3) where ∂S ≡ (d/dt)S and ∂S' ≡ (d/dt)S', ≡ ∂SL = ∂S'L + ω x L (2.1) so Newton (I.2.1) says (∂S'L) + ω x L = N . (I.2.2) Note that (∂S'L) is not a "natural" object in the sense of Section 1.8 since L is a Frame S object, but the time derivative is taken in Frame S', the body frame. Now consider, (∂S'L) = (∂S'[Iω]) = (∂S'I)ω + I(∂S'ω) . (I.2.3) We noted above that the inertia tensor is constant in Frame S' so we expect that (∂S'I) = 0. To make sure, we examine the components of (∂S'I) in Frame S', [(∂S'I)]'ij = ∂S'(I)'ij // commutation theorem (1.11.1) applied to a tensor = ∂t(I')ij // (1.10.1) = 0 . // shown below (I.1.6) (I.2.4) If [(∂S'I)]'ij = 0 for all components, then (∂S'I) = 0 as a tensor statement. We also know from (2.6) that ∂S'ω = ∂Sω = ∂tω = . (2.6) Therefore (I.2.3) reads, (∂S'L) = (∂S'[I ω]) = (∂S'I) ω + I (∂S'ω) = 0 + I = I . (I.2.5) Inserting this into (I.2.2) then produces our vector equation of motion for ω, N = I + ω x (Iω) . (I.2.6) Taking components in Frame S' (the body frame) one finds (N)'i = (I)'ij ()'j + εijk(ω)'j(I)'ka (ω)'a . (I.2.7) In what follows we shall always solve (I.2.6) in Frame S' components because only in Frame S' can the inertia tensor be diagonalized in a manner that applies for all time t. I.3 Diagonalization of the Inertia Tensor in Frame S' So far we have not paid much attention to the manner in which the Frame S' basis vectors are selected, and so far we have called them e'i. One is of course free to define new orthonormal Frame S' basis vectors e"i according to e"i ≡ Qij e'j or | e"i> ≡ Qij | e'j> <e'k | e"i> = Qik where Q is some arbitrary rotation matrix (Q-1 = QT). Then in this new basis the inertia tensor is (I)"ij = <e"i| I | e"j> = <e"i |e'a><e'a| I |e'b><e'b| e"j> = = Qia(I)'abQjb = [Q I Q-1]ij // (I)" = Q I Q-1 Note that (I)'ij = Σαmα [ r'α2δij - (rα)'i(rα)'j ] is a symmetric matrix, (I)'ij = (I)'ji. A staple matrix theorem is that any real symmetric matrix M can be brought to diagonal form by a "similarity transformation" with some real orthogonal matrix S, so Λ = S-1MS where Λ is diagonal. There is a specific method for determining a suitable matrix S that does this diagonalization task. If we then select Q = S-1, the matrix (I)" = Q I Q-1 will be diagonal in the e"i basis. The basis vectors e"i are called "the principal axes" of the rigid body in question, and Q = S-1 is "the principal axis transformation". Note that the principle axes are fixed to a rigid body and do not change relative to the body as the body is reoriented (like the center of mass, and unlike the center of gravity, as both discussed in Appendix D). Realizing this fact, we might as well start off selecting e'i to be a set of principle axes. From now on, then, we shall assume that the e'i have been so selected and thus the inertia tensor I is diagonal in Frame S' components. (Remember that the Frame S inertia tensor I' is identically zero.) Comments: In general, any real symmetric matrix M can be diagonalized by some real orthogonal matrix S to become the diagonal matrix Λ = S-1MS where then Λ = diag(λ1,λ2,λ3). Suppose the matrix M has some eigenvectors vi with eigenvalues μi so that Mvi = μivi. Then since M = S-1MS one gets, (SΛS-1)vi = μivi Λ(S-1vi) = μi(S-1vi) Λwi = μiwi . // wi ≡ S-1vi But the eigenvalues of a diagonal matrix Λ are those diagonal elements, so we identify λi = μi and conclude that when M is diagonalized to Λ, the diagonal elements are the eigenvalues of M. The eigenvalues can be found by writing the eigenvalue equation as (M - μi1)vi = 0. If det(M-μi1) ≠ 0 one can invert this equation to conclude that vi = 0 which is wrong, so one must have det(M-μi1) = 0. This so-called "secular equation" or "characteristic equation" is then a cubic in variable μi which can be solved for three values of μi (the three roots of the secular equation). In the complex world, the claim is that any Hermitian matrix H can be diagonalized by a unitary transformation U. A matrix H is Hermitian if H†= H( HT* = H), and a matrix U is unitary if U†U = 1 so U† = U-1. Then one has Hd = U-1HU. As above, the diagonal elements of Hd are the eigenvalues of H. It is easy to show that the eigenvalues of H (which are the diagonal elements of Hd) are real: Hd* = U*-1H*U* = UT* HT*U*-1T = U† H†U-1† = U-1HU = Hd. In quantum mechanics, all observable quantities X are represented by Hermitian operators (matrices) and the diagonal elements of the diagonalized Xd (the eigenvalues of X) are the possible physical values the observable can take. Given that the e'i are selected as a set of principal axes for a rigid body, the diagonal elements of the diagonal inertia tensor can be written I'i and are called the principal moments of inertia. So (I)'ij = I'i δij . (I.3.1) We maintain a prime on I'i as a reminder that these are associated with Frame S' where the inertia tensor I has been made diagonal by selection of the e'i basis vectors of Frame S'. I.4 Rigid Body Equations of Motion : Part 2 We now insert the diagonal form (I.3.1) for (I)'ij into the equations of motion (I.2.7) to get (N)'i = I'i δij ()'j + εijk(ω)'j I'k δka (ω)'a = I'i ()'i + εijk(ω)'j I'k(ω)'k = I'i ()'i + εijk(ω)'j(ω)'k I'k . (I.4.1) Setting i =1 gives (N)'1 = I'1 ()'1 + ε1jk(ω)'j(ω)'k I'k = I'1 ()'1 + ε123(ω)'2(ω)'3 I'3 + ε132(ω)'3(ω)'2 I'2 = I'1 ()'1 + (ω)'2(ω)'3 I'3 - (ω)'3(ω)'2 I'2 = I'1 ()'1 - (ω)'2(ω)'3 (I'2 - I'3) . (I.4.2) Examining the other components (or just doing cyclic permutations), we may state the three component rigid-body equations of motion as follows (all components are Frame S' components), (N)'1 = I'1 ()'1 - (ω)'2(ω)'3 (I'2 - I'3) (N)'2 = I'2 ()'2 - (ω)'3(ω)'1 (I'3 - I'1) (N)'3 = I'3 ()'3 - (ω)'1(ω)'2 (I'1 - I'2) . (I.4.3) This is a non-linear system of first-order ODE's and they appear in Goldstein p 158 (5-34) or GPS p 200 (5.39'). At this point Goldstein is using our swap notation so there are no primes on any of the quantities. For torque-free motion one then has I'1 ()'1 = (ω)'2(ω)'3 (I'2 - I'3) I'2 ()'2 = (ω)'3(ω)'1 (I'3 - I'1) I'3 ()'3 = (ω)'1(ω)'2 (I'1 - I'2) . (I.4.4) Goldstein p 159 (GPS p 201) notes that equations (I.4.4) can in principle be analytically solved in terms of elliptic functions, but that the resulting complicated solutions are not very enlightening. We have seen what such solutions look like for the spherical pendulum in (C.5.13), and the elliptic forms appear again in (I.7.17) for the treatment of the gravity top. I.5 Zero-torque motion of a Rigid Body : Ellipsoids and Poinsot One might naively think that an arbitrary rigid body floating in space which is initialized to rotate about some rotation axis with some ω might continue to do so. If it happens that the initialized ω points in the direction of the largest or the smallest principal axis of inertia of the rigid body, this is in fact what happens. Otherwise the rigid body tumbles in a complicated manner in order to satisfy equations (I.4.4). Suppose the principle moments of inertia are such that I'3 > I'2 > I'1 . If rotation is initialized such that ω = ωe'2 so that only (ω)'2 ≠ 0, one finds that the three equations (I.4.4) are in fact satisfied by solution (ω)'2 = constant. However, as shown for example in Marion p 399 (T&M p 460), the tiniest perturbation causes the rigid body to lose its simple motion and to tumble in a complicated manner. This is traditionally demonstrated by initializing a rubber-banded book or tennis racket to rotate about its middle principle axis while tossing the object a few feet into the air. Even in this short time, the I'2 axis rotation cannot be maintained and the object tumbles. In contrast, rotation is stable about the other two principal axes. In lieu of the complicated elliptic solutions noted above, there are various geometric constructions involving ellipsoids that are used to interpret torque-free motion of a rigid body. If the e'i are the principal axes of a rigid body in Frame S', one may write (I.1.9) and (I.1.4) squared as follows, 2T = ωIω = (I)'i(ω)'i2 1 = + + // inertia ellipsoid L2 = Iω Iω = (I)'i2(ω)'i2 1 = + + . // L-ellipsoid (I.5.1) A torque-free isolated system will have fixed values for rotational kinetic energy T and angular momentum L2. One can then regard each of the above equations as an axis-aligned ellipsoid centered at the origin of ω-space with axes e'i. The denominators shown are the squares of the three semi-axes of each ellipsoid. The rigid body of interest is aligned with the inertia ellipsoid since it is aligned with the e'i basis vectors in Frame S'. If I'3 is the smallest principal moment of inertia (suggesting prolateness of the object in that direction), then the longest axis of the inertia ellipsoid lies in the (ω)'3 direction. The inertia ellipsoid is sometimes called Poinsot's ellipsoid. Goldstein defines ρ ≡ (1/)ω as a rescaled version of ω and then one has 1 = ρIρ = (I)'i(ρ)'i2 1 = (I)'1(ρ')12 + (I)'2(ρ')22 + (I)'3(ρ')32 . This is the "official" inertia ellipsoid, but we shall not do this rescaling and shall refer to the first line of (I.5.1) as the inertia ellipsoid. Since rigid body problems do have solutions, we know that the two ellipsoids of (I.5.1) must intersect. In general, the intersection of two axis-aligned ellipsoids is a set of two non-planar curves which are mirror images of each other. In the solution of a torque-free rigid body problem the tip of the ω vector must move along one of these intersection curves. One can thus regard such a curve as a closed "orbit" for the ω vector in Frame S'. This is somewhat analogous to finding the orbit r(θ) for the polar coordinate position r of a planet without actually solving for the detailed r(t) and θ(t) of the planet's motion. Finding the equation of a curve that is the intersection of two ellipsoids is not easy, even when the ellipsoids are axis-aligned, but it is easy to graphically display the curved intersection orbits. In the Maple code below we enter three different ellipsoids Red, Blue and Green, where the Green one is a sphere. One should think of [x,y,z] as [(ω)'1 ,(ω)'2,(ω)'3] = ω in Frame S' components. We then plot the Red and Green ellipsoids, (I.5.2) One can see that the intersection curves are distinctly non-planar and warped like a potato chip. Nevertheless, in this example the length of the ω vector is constant on each orbit curve, since one ellipsoid is a sphere. (The Green sphere could be the inertia ellipsoid of a tumbling cube.) We now look at the Red and Blue ellipsoids, (I.5.3) Again, the intersection curves are non-planar. We then ask: it is possible that ω has a constant magnitude on these Red/Blue intersection curves? If so, then the curves must also lie on a sphere. We now add the Green sphere into the above picture : (I.5.4) Although one knows that the Red/Green intersection curves are warped, one can see that their shape is never going to match the shape of the Red/Blue intersection curves no matter what radius is selected for the Green sphere: We therefore reach a conclusion that is perhaps obvious: Fact: When equations (I.4.4) are solved, the orbit of the solution vector ω(t) is such that the magnitude ω is (in general) not constant. Magnitude ω is constant only if the tip of ω(t) does conical motion, as vector a does in (1.6.2) where the tip goes in a circle. But the warped orbits shown above in general are not circles. In our final example, suppose the rigid body is axisymmetric so I'1 = I'2. We adjust the example above so that (I.5.5) Here we selected the green sphere radius so the Green sphere is just visible. It seems clear in this situation that the Red/Blue intersection curves are no longer warped, but are just circles perpendicular to the z axis on which ω moves, this z axis being the (ω)'3 axis. As we shall see below in Fig (I.6.14), in this case the ω vector does do conical motion where its tail lies at the center of the ellipsoids (the origin of ω space) and its tip goes around one of the Blue/Red intersection circles of (I.5.5). The ω vector moves at a uniform rate as shown below, and this motion of ω is called precession. An ellipsoid intersection curve of the type seen in the above examples is called a polhode (a pole path), polhode: mod. f. Gr.  pole +  way, path (Poinsot 1852) // OED2 For a given choice of the (I)'i moments, if one keeps T fixed and varies L in (I.5.1), one generates a family of polhode intersection curves on the inertia ellipsoid, as shown here (Arnold and Maunder p 111, moments of inertia are called A,B,C ) (I.5.6) As already noted, in ω-space, the solution vector ω(t) has its tail glued to the center of the inertia ellipsoid while its tip rides along a polhode. The drawing shows that for the largest and smallest inertia axes, the extremal polhodes are circular, suggesting a possible slightly circular motion of a rigid body axis set to rotate close to one of these axes. There are no such circles for the intermediate axis, which is why that axis is unstable against perturbation as was noted above. Any tiny variation away from the exact axis point causes the solution to run off on a large polhode which encircles the ellipsoid. The conclusion so far is that, for torque-free motion of a rigid body, the ω(t) vector in Frame S' (the body frame) moves on some relatively simple closed polhode path which is in general non-planar. What does one of these potato-chip warped polhode paths look like when viewed from inertial Frame S? Perhaps surprisingly, one thing we do know is that in Frame S the path of the ω vector tip lies in a plane, so the potato chip polhodes are de-warped in their transformation from Frame S' to Frame S. The reason for this planarity is simple. In Frame S we know that the vector L is fixed in space (since no torque), and one has ω = (1/L) L ω = (1/L) (Iω) ω = (1/L) ω I ω = (1/L)(2T) = 2T/L . (I.5.7) Now in E3 space the equation r = d describes a plane with normal vector lying distance d from the origin (d is the closest distance from plane to origin). Thus, in ω-space the vector ω always lies on a plane with normal which lies distance 2T/L from the origin. Since this plane is fixed in Frame S, it is called "the invariable plane". At any time t, the solution ω(t) lies both on the inertia ellipsoid and on the invariable plane, so the inertia ellipsoid viewed in Frame S must be tangent to the invariable plane. Here is a picture where we have drawn pointing down, (I.5.8) If the ellipsoid were an ellipse and if the above were a 2D picture, the ellipse would be unable to roll on the "invariable line", since any such rolling would change the height of the ellipse center above the plane. But in the 3D picture with an ellipsoid, such rolling is possible. If the ellipsoid were axisymmetric, it would roll around the vertical L axis, thus inscribing a circle on the invariable plane. That, we shall see, is the solution found in the next section and shown in (I.6.30). If the ellipsoid is more general in shape, then the rolling motion of the ellipsoid produces some complicated path going round and round on the invariable plane. This path is called a herpolhode (a snaky pole path) and is generally a path which does not close on itself, but which is confined between two radii in the invariable plane. Here is a sample numerically-generated herpolhode (Arnold and Maunder p 111 ) (I.5.9) If the ellipsoid in Fig (I.5.8) were covered with ink and if the invariable plane were paper, the rolling ellipsoid would inscribe the herpolhode on the invariable plane. Conversely, if the ellipsoid were paper and the invariable plane were coated with ink, the invariable plane would inscribe on the ellipsoid one of the closed, non-planar polhode paths. In general, this geometric view of the ω vector motion is called "Poinsot's construction" (Louis Poinsot 1777 –1859). It is briefly described in Goldstein p 159 (GPS p 201). Here is an excellent ellipsoid animation showing both polhode and herpolhode for a case where the latter is also a closed curve: https://www.youtube.com/watch?v=BwYFT3T5uIw . I.6 Zero-torque motion of an Axisymmetric Rigid Body (Rigid Rotor) Finding the Frame S' components of ω Here we specialize to rigid objects for which two of the principal moments of inertia are equal, and we shall take the two equal moments to be I'1 = I'2. A rigid object which is a solid of revolution falls into this class, where the symmetry axis will be the axis with I'3. The other two orthonormal principal axes can be selected in any manner to lie in the plane perpendicular to the symmetry axis. For a pancake, the symmetry axis has the largest principal moment, while for a pencil it has the smallest principal moment. Less symmetric objects can have I'1 = I'2 as well, such as a rod whose cross section is a regular polygon with an even number of faces. For such a rod, the two transverse principal axes can be any two orthonormal axes which are also symmetry axes of the rod cross section. All objects having the same values of I'3 and I'1 = I'2 tumble in exactly the same manner, regardless of their shape. Example 1: Square centered at the origin of uniform density ρ : (I.6.1) We know that Ixy = - ρ∫dxdy (xy). In the left picture one can see that the positive xy contributions from quadrants I and III exactly offset the negative xy contributions from quadrants II and IV. Thus Ixy = 0 and the axes shown are principal axes. The exact same is true for the other two figures because the shape in each quadrant is the same. So any perpendicular pair of axes through the enter of a square can serve as principal axes. Similarly for a cube, any three orthogonal axes through the center point can be principal axes. Because the inertia tensor has only one set of eigenvalues (which are the principal moments), those moments are the same for any legal set of principal axes. Example 2: Hexagon centered at the origin of uniform density ρ : (I.6.2) In this case the previous example's argument holds for the left two pictures, so those axes can be taken to be principal axes. For the orientation on the right, the QI and QII shapes are not the same. In this case only the left two pictures show principal axes. If e'3 is the symmetry axis and I'1 = I'2, the last of equations (I.4.4) reads 0 = I'3()'3 , so (ω)'3 = constant. The other two equations in (I.4.4) are then easily solved. Writing them out, ()'1 = (ω)'3 (ω)'2(I'1 - I'3)/ I'1 = - Ω'(ω)'2 where Ω' ≡ (ω)'3 (I'3-I'1)/I'1 ()'2 = (ω)'3 (ω)'1(I'3 - I'1)/I'1 = Ω'(ω)'1 . (I.6.3) Then ()'1 = - Ω'()'2 = - Ω'2 (ω)'1 ()'1 + Ω'2 (ω)'1 = 0 . (I.6.4) A simple solution to his harmonic motion ODE is as follows, (ω)'1 = A' cos(Ω't) (ω)'2 = - ()'1 / Ω' = A' sin(Ω't) . (ω)'3 = K (I.6.5) Without loss of generality, we take A' > 0 and so A' = . This solution describes "conical motion" (as for vector a in (1.6.2)) of the ω vector where the cone half-angle α is determined by sinα = A'/ω cosα = K/ω tanα = A'/K . (I.6.6) For (ω)'3 = K > 0, α lies in the range (0,π/2), and for (ω)'3 = K < 0, α lies in the range (π/2,π) . Thus we have this simple precession solution for the Frame S' components of the ω vector for zero-torque rotation of an axisymmetric rigid body : (ω)'1 = ωsinα cos(Ω't) Ω' = ωcosα (I'3-I'1)/I'1 (ω)'2 = ωsinα sin(Ω't) (ω)'3 = ωcosα . α = cone half-angle, 0 ≤ α ≤ π (I.6.7) The direction of the precession depends on the sign of cosα and on the sign of I'3-I'1. For a "fat" rotor (oblate), the symmetry axis will have the larger moment, so I'3 > I'1. For α < π/2 this means Ω' > 0, and the precession is CCW about the cone axis as viewed from the circle end of the cone. One can then ask about the behavior of the L vector in Frame S'. Since (L)'i = I'i(ω)'i one has, (L)'1 = I'1ωsinα cos(Ω't) Ω' = ωcosα (I'3-I'1)/I'1 (L)'2 = I'1ωsinα sin(Ω't) (L)'3 = I'3ωcosα . (I.6.8) Notice that for a "very fat rotor" (oblate) where I'3 >> I'1, most of the angular momentum is in the (L)'3 component. But for a "very thin rotor" (prolate) where I'3 << I'1 , most of the angular momentum is in the transverse precession components (L)'1 and (L)'2 . From this last set of equations one finds that L2 = ω2 ( I'12sin2α + I'32cos2α ) (I.6.9) where L is the magnitude of L. The kinetic energy T of the rigid body is determined from (I.6.8) to be, 2T = ω I ω = I'i(ω)'i2 = I'1[ ωsinα cos(Ω't)]2 +I'2[ ωsinα sin(Ω't)]2 + I'3 [ ωcosα]2 = I'1 ω2sin2α + I'3ω2cos2α = ω2 ( I'1sin2α + I'3cos2α ) . (I.6.10) To summarize, L2 = ω2 ( I'12sin2α + I'32cos2α ) T = (1/2) ω2 ( I'1sin2α + I'3cos2α ) . (I.6.11) If ω and α are specified, then L and T are determined by these equations. Conversely, if L and T are specified one can determine values for α and ω. Reader Exercise: Show that the solutions for α and ω are given by tan2α = (I'3/I'1) ω2 = + (I.6.12) Returning now to the L equations (I.6.8), (L)'1 = I'1ωsinα cos(Ω't) Ω' = ωcosα (I'3-I'1)/I'1 < 0 (I.6.8) (L)'2 = I'1ωsinα sin(Ω't) (L)'3 = I'3ωcosα . one sees that the L vector (in Frame S') precesses on a cone at the same rate Ω' that ω precesses in (I.6.7). However, for L the cone half-angle β is determined by sinβ = I'1ωsinα /L cosβ = I'3ωcosα/L tanβ = (I'1/I'3) tan α . (I.6.13) In all the figures bellow we assume that 0 < α < π/2 so cosα > 0. The size of β relative to α depends on whether the inertia moment ratio (I'1/I'3) is > 1 or < 1. If the rotor is "thin" (prolate), then I'1 > I'3 and (I'1/I'3) > 1. Therefore β > α and Ω' = ωcosα (I'3-I'1)/I'1 < 0. This is shown in the left figure below. If the rotor is "fat" (oblate), then I'1 < I'3 and (I'1/I'3) < 1. Therefore β < α and Ω' = ωcosα (I'3-I'1)/I'1 > 0. This is shown in the right figure below. (I.6.14) In either case, the ω and L vectors rotate in the same direction. The dizzy Frame S' Observer of course sees the rigid body rotor at complete rest with its symmetry axis ("figure axis") pointing in the ' = e'3 direction ("up"). Frame S and any objects at rest in Frame S are seen to be violently rotating about this Frame S' Observer. The vector L is the Frame S angular momentum of the rigid body, but we are just showing it in Frame S' coordinates. The Frame S' angular momentum L' of the rigid body is of course 0. The Observer can note that -ω is the angular velocity of Frame S relative to Frame S'. The Figures are not particularly interesting, but they do show the solution to the problem in Frame S' coordinates. Soon we shall show that β = θ, an Euler angle of the rigid body in Frame S. The behavior of the Rigid Rotor in Frame S : solving for the Euler angles Of much more interest is what things look like to an Observer in inertial Frame S. This takes a bit more work. First, we shall assume that the position of the symmetric rigid rotor is described by the Euler angles (φ,θ,ψ) of Fig (H.1.5), so in Frame S one has (I.6.15) Frame S has axes x,y,z while Frame S' has axes x',y',z'. Recall from (H.3.12) the transformation for a vector r going from Frame S to Frame S', x' = (cosψcosφ - sinψcosθsinφ) x + (cosψsinφ + sinψcosθcosφ) y + sinψsinθ z y' = (- sinψcosφ - cosψcosθsinφ) x + (-sinψsinφ + cosψcosθcosφ) y + cosψsinθ z z' = sinθsinφ x - sinθcosφ y + cosθ z . (H.3.12) In Frame S where N = and N = 0, L must be a constant vector. We arbitrarily put this in the +z direction, so that L = L in Frame S. Applying the above transformation to L instead of r gives, (L)'1 = sinθsinψ L (L)'2 = sinθcosψ L (L)'3 = cosθ L . (I.6.16) But our Frame S' problem solution (I.6.8) with (I.6.13) gives (L)'1 = Lsinβ cos(Ω't) Ω' = ωcosα (I'3-I'1)/I'1 (L)'2 = Lsinβ sin(Ω't) (L)'3 = Lcosβ . (I.6.17) Comparison of these last two equation sets requires that, sinθsinψ = sinβ cos(Ω't) so: cos(π/2-ψ) = sinψ = cos(Ω't) sinθcosψ = sinβ sin(Ω't) sin(π/2-ψ) = cosψ = sin(Ω't) Ω't = π/2-ψ cosθ = cosβ . (I.6.18) The solution to these equations is θ = β ψ = - Ω't +π/2 = - Ω' . (I.6.19) This says that Euler angle θ is a constant, and also, from (I.6.13), sinθ = I'1ωsinα /L cosθ = I'3ωcosα/L tanθ = (I'1/I'3) tan α . (I.6.20) We now rewrite (I.6.7) as (ω)'1 = ωsinα sinψ (ω)'2 = ωsinα cosψ (ω)'3 = ωcosα . (I.6.21) Now if we set = 0 in our Frame S' Euler-angle evaluation of the (ω)'i given in (H.4.7), we find (ω)'1 = sinθsinψ // ψ = Ω't, = Ω' (ω)'2 = sinθcosψ (ω)'3 = cosθ + = cosθ - Ω' . (H.4.7) (I.6.22) Comparison of these last two equation sets requires that, sinθ = ωsinα cosθ - Ω' = ωcosα . (I.6.23) From the first line above and then using (I.6.20) one finds, = ωsinα/sinθ = L/I'1 φ(t) = (L/I'1) t . (I.6.24) We have now solved for the behavior of all the Euler angles which describe the orientation of our rigid rotor in Frame S: θ = tan-1[ (I'1/I'3) tanα ] // θ(t) = constant ψ(t) = - Ω't + π/2 = - Ω' // Ω' = ωcosα (I'3-I'1)/I'1 φ(t) = (L/I'1)t = L/I'1 . (I.6.25) With regard to Fig (I.6.15), the rigid rotor spins at rate -Ω' about its symmetry z' axis. While this happens, the symmetry axis precesses CCW at rate = (L/I'1) about the z axis. There is no nutation since the polar angle θ is a constant. The precession is CCW because we assumed that L = L where L = |L| . (If the rotor is a cube, since I'3 = I'1 one has = 0 since Ω' = 0. ) The behavior of the Rigid Rotor in Frame S : solving for ω Recall the Frame S' solution for ω shown in (I.6.21), (ω)' = . // sinψ = cos(Ω't) and cosψ = sin(Ω't) (I.6.21) To obtain the Frame S components of ω we have Maple compute ω = R-1(ω)' : where W = ω . The resulting vector can be simplified by combining the trig terms: The final result is then ω1 = ωsin(θ-α)sinφ ω2 = - ωsin(θ-α)cosφ ω3 = ωcos(θ-α) . (I.6.26) Reader Exercise. (1) Use Ω' = ωcosα (I'3-I'1)/I'1 from (I.6.7) and tanθ = (I'1/I'3) tan α from (I.6.20) to show that Ω' = ω sin(α-θ)/sinθ. (2) Using (I.6.25), show that (I.6.26) for ωi agrees with (H.4.12) for ωi. (I.6.27) With φ = (L/I'1)t from (I.6.24), (I.6.26) can be written, ω1 = ω sin(θ-α) sin[(L/I'1)t ] ω2 = - ω sin(θ-α) cos[(L/I'1)t ] ω3 = ω cos(θ-α) . (I.6.28) or ω1 = ω sin(θ-α) cos[(L/I'1)t -π/2] ω2 = ω sin(θ-α) sin[(L/I'1)t -π/2] ω3 = ω cos(θ-α) . (I.6.29) If 0 < θ-α < π, then sin(θ-α) > 0. In this case, equations (I.6.29) describe (in Frame S) the CCW precession at rate(L/I'1) of the vector ω. Here for the prolate case is the famous picture traditionally used to torture students of rotational dynamics, where one should momentarily ignore the red cone : (I.6.30) As just noted, ω does CCW conical motion at rate (L/I'1) on a cone of half-angle θ-α, shown in blue. The rigid rotor's symmetry axis meanwhile moves on the black θ cone and precesses along with ω at the same rate (L/I'1). The is really the end of the Frame S story since it tells how the rotor figure axis moves, how its ω vector moves, and how L is fixed in the vertical direction since this is a torque-free system in Frame S. One can (if one wants) superpose the red body cone from our Frame S' picture (I.6.14). Since there is only one ω vector, the body cone must be at a position like that shown. As blue ω rotates around on its fixed space cone, the red body cone must roll without slipping around the perimeter of the blue cone to maintain ω on the intersection of the two cones. Notice that the red direction arrows on the body cone match the rotation sense of the red arrow in the left drawing in (I.6.14). To see why this corresponds to rolling without slipping, one must physically play with two cones. A substitute is two coins where one holds the left one fixed and rotates the right one about the boundary of the left one. Only then does one understand the red arrows! It is too difficult to put into words. One can see that when the red cone is made narrower (smaller α), the angular rate of motion of ω on the red cone perimeter increases relative to the angular rate of ω on the blue cone. The blue rate is (L/I'1) while the red rate is Ω' = ωcosα (I'3-I'1)/I'1. At this point, as an illustration of the above figure, the reader is invited to view Eric Johnson's short animation https://www.youtube.com/watch?v=s9wiRjUKctU showing the motion of a prolate rotor. One sees how his purple L vector (called H) says fixed, and how the blue ω vector moves on its cone, and how the figure axis moves on a larger cone. An American football in flight is a good example of a prolate rotor, see Horn and Fearn. On the other hand, if the rotor is fat (oblate) in shape with I'1< I'3, then Ω' > 0 and θ < α. The corresponding figure is shown below, where the blue and black cones of (I.6.30) have not changed : (I.6.31) Now the red arrows on the body cone have the same directional sense as the right drawing in (I.6.14) . Reader exercise: Use two coins or two drink coasters to verify the direction of the red arrows. Eric Johnson's oblate rotor animation is here: https://www.youtube.com/watch?v=PDLXVSkDFVk . Again the purple L = H vector is fixed, but now the figure axis cone lies outside the ω cone. See https://www.youtube.com/watch?v=WUkUL3Hp67A for an animation of the prolate case body cone rolling around the space cone. Reader Question: Since L = Iω, and since L = constant, one can write ω = I-1L . Why is ω not constant? Answer: In Frame S, the inertia tensor I (relative to the fixed Frame S axes) constantly changes with time because the integration in (I.1.3) is a function of time, as if the rotating rotor were an object constantly changing shape. Another answer: I = R-1(I)' R and R = Rz(-ψ(t)) Rx(-θ(t)) Rz(-φ(t)) varies with time. In fact one has for the axisymmetric case (Reader Exercise), I(t) = Rz(φ)Rx(θ)Rz(ψ)(I)'Rz(-ψ(t)) Rx(-θ(t) Rz(-φ(t)) = Rz(φ)Rx(θ)(I)'Rx(-θ(t) Rz(-φ(t)) . The Earth as an oblate rigid rotor An ideal model of the Earth has it being a rigid rotor that is slightly oblate due to the centrifugal force on the slightly elastic matter of the Earth (and its water) pulling matter out to a radius larger than the average radius of the Earth as one approaches the equator. As Goldstein notes on p 163 (GPS p 207), the moment ratio appearing in Ω' is about .0033 which predicts a precession rate of Ω' ≈ (ω)'3/300 or T ≈ 300 days. Various inadequacies in this idealized model (Earth is not rigid, nor is it an exact oblate spheroid) are blamed for the fact that the precession period is more like 434 days. Due to this precession the North pole moves in a circle of radius about 6 m. The precession cone half angle is a tiny 0.2 arc seconds, meaning 0.2/3600 of one degree. This effect is called the "free precession of the Earth" and is also known as the Chandler wobble (S.C. Chandler 1891). Sometimes this precession is called a nutation since it is superposed on longer timescale precessions due to the torque of the Sun and Moon acting on the Earth as described below. Note that the free precession effect of the spinning Earth would occur if the Earth were completely isolated in space. It has nothing to do with external torques acting on the Earth. In fine detail, there are many factors which contribute to the Earth's axis wobble which in its aggregate is called the Chandler wobble. For an ideal axisymmetric oblate Earth, the polhode path should be circular as in Fig (I.5.5) or (I.6.31). Here are some real world polhodes (Zemstov), where each small square is 3 meters on a side (it seems that the squares should be 0.1" rather than 0.01"), (I.6.32) I.7 Rigid Body with External Torque: Spinning Symmetric Top In this problem Frames S and S' have their origins co-sited at the top's point of contact with a "table". Frame S is the inertial frame of the table, with z pointing up. Frame S' is embedded in the spinning top. We assume that the top is symmetric and is spinning about its symmetry axis e'3 = ' and then I'1 = I'2. If we put a red paint dot on the top at some value of ψ, we can describe the position of the top at some time t by the Euler angles φ,θ,ψ of Section H.1 : (I.7.1) Gravity g = -g creates a downward force F = Mg = -Mg acting on the center of mass of the top, which lies at point rcms on the top's symmetry z' axis, some distance rcms up from the pivot point. Since F is parallel to while rcms is parallel to ', the torque N = rcms x F (referred to the common origin) must be perpendicular to both and '. The locus of points perpendicular to and ' is the line of nodes in the figure, the intersection of the two discs, so N points toward the viewer along the ' axis. The fact that the torque is perpendicular to the and ' axis reveals that the angular momenta associated with angles φ and ψ are constant. Those angular momenta are Lz = L3 for φ, and (L)'z = (L)'3 for ψ. In Lagrangian language, this means that the coordinates φ and ψ do not appear in L and are therefore cyclic, so their generalized momenta pψ = (L)'z and pφ = (L)z are constants of the motion. Following the Lagrangian method, note first from (I.1.9) the expression for kinetic energy T in Frame S' components, T = (1/2) (ω)'i(ω)'j(I)'ij = (1/2)(ω)'i2 I'i = (1/2)[ (ω)'12 + (ω)'22 ] I'1 + (1/2)(ω)'32 I'3 . (I.7.2) Recall from (H.4.7) that (ω)'1 = sinθsinψ + cosψ (ω)'2 = sinθcosψ - sinψ (ω)'3 = cosθ + // Frame S' (H.4.7) (I.7.3) so that (ω)'12 + (ω)'22 = [ sinθsinψ + cosψ]2 + [ sinθcosψ - sinψ]2 = 2 sin2θ + 2 . // cross terms cancel (I.7.4) Thus T = (1/2)(2 sin2θ + 2) I'1 + (1/2)( cosθ + )2 I'3 . (I.7.5) If the zero of potential is set at z = 0, the Lagrangian L is then L = T - V = (1/2)(2 sin2θ + 2) I'1 + (1/2)( cosθ + )2 I'3 - Mgrcmscosθ (I.7.6) confirming that φ and ψ are cyclic. The two constant generalized momenta are then (they are constant since ∂L/∂φ = 0 and ∂L/∂ψ = 0), pψ = ∂L/∂ = ( cosθ + )I'3 ≡ aI'1 // = (ω)'3I'3 (I.7.7) pφ= ∂L/∂ = sin2θ I'1 + ( cosθ + )cosθ I'3 = (sin2θ I'1 + cos2θ I'3) + cosθ I'3 ≡ bI'1 , where constants a and b are defined as shown. Since (ω)'3 = cosθ + , one sees from the pψ line above that, (ω)'3 = (aI'1/I'3) so (ω)'3 is a constant, consistent with the claim above that there is no torque about the z' axis. Exercise: Show directly (without using the Lagrangian) that the constants of the motion (L)'z and (L)z are the same as the generalized momentum expressions pψ and pφ stated above. (a) From L = Iω we know that (L)'z = (L)'3 = I'3 (ω')3 = I'3 ( cosθ + ), thus (L)'z = pψ . QED1 (b) Write (L)z = L = L [(sinψsinθ)' + (cosψsinθ)' + (cosθ) ' ] using (H.3.18). Then (L)z = sinψsinθ(L')1 + cosψsinθ(L')2 + cosθ(L')3 = sinψsinθ I'1(ω)'1 + cosψsinθ I'2(ω)'2 + cosθ I'3(ω)'3 = sinψsinθ I'1[ sinθsinψ + cosψ] + cosψsinθ I'2[ sinθcosψ - sinψ] + cosθ I'3[ cosθ + ] = [sin2θI'1 + cos2θI'3 ] + cosθ I'3 since I'1 = I'2, thus (L)z = pφ . QED2 We are closely following Goldstein page 165 (GPS p 211) where the constants of the motion (L)'z and (L)z are replaced by constants a and b. The equations of interest (I.7.7) are then I'3 cosθ + I'3 = aI'1 (sin2θ I'1 + cos2θ I'3) + cosθ I'3 ≡ bI'1 . (I.7.8) Multiply the first equation by cosθ and subtract 2nd - 1st to get (sin2θ I'1 + cos2θ I'3) - I'3 cos2θ = bI'1 - cosθ aI'1 sin2θ I'1 = bI'1 - cosθ aI'1 sin2θ = b - cosθ a so = . (I.7.9) Now put this into the first equation of (I.7.8) to get I'3[( b - cosθ a)/sin2θ] cosθ + I'3 = aI'1 ( b - cosθ a)/sin2θ] cosθ + = a (I'1/I'3) = a (I'1/I'3) - cosθ( b - cosθ a)/sin2θ . (I.7.10) Once θ(t) is determined as outlined below, one can integrate (I.7.9) and (I.7.10) to get φ(t) and ψ(t). A third constant of the motion is the total energy of the top, where T was stated in (I.7.5), E = T + V = (1/2)(2 sin2θ + 2) I'1 + (1/2)( cosθ + )2 I'3 + Mgrcmscosθ . (I.7.11) The second term in E is just a constant from (I.7.7), (1/2)( cosθ + )2 I'3 = (1/2)(ω'3)2 I'3 = (1/2)(aI'1/I'3)2 I'3 = (1/2)a2(I'12/I'3) = K (I.7.12) so we define E' = E - K as a re-zeroed energy to get, E' = (1/2)(2 sin2θ + 2) I'1+ Mgrcmscosθ . (I.7.13) Multiply by (2/I'1), (2E'/I'1) = (2 sin2θ + 2) + (2Mgrcms/I'1) cosθ or α = (2 sin2θ + 2) +β cosθ where α ≡ (2E'/I'1) β ≡ (2Mgrcms/I'1) or α sin2θ = 2 sin4θ + sin2θ 2 + βcosθsin2θ . (I.7.14) From (I.7.9) replace 2 sin4θ by (b-acosθ)2, α sin2θ = (b-acosθ)2 + sin2θ 2 + βcosθsin2θ and rearrange to get sin2θ 2 = sin2θ(α - βcosθ) - (b-acosθ)2 . (I.7.15) Set u = cosθ so = -sinθ and then sin2θ 2 = 2, so the ODE becomes 2 = (1-u2)(α-βu) - (b-au)2 = α - u2α - βu +βu3- b2 +2abu - a2u2 = βu3 - (α+a2)u2 + (2ab-β)u + (α-b2) = β(u-A)(u-B)(u-C) (I.7.16) where one can write A,B,C in terms of a,b,α,β. This is a solvable non-linear first-order ODE. Write du/dt = and so dt = du/ and then t(u) = (1/) !Syntax Error, Idx . (I.7.17) This is the exact same elliptic integral we encountered with the spherical pendulum in (C.5.13). One does the integral to get t(u) = (1/) [ f(A,B,C,u) - f(A,B,C,u0) ] and then one "inverts" to get u = u(t) and thus one has found θ(t) = cos-1u(t). Functions φ(t) and ψ(t) are then found by integrating (I.7.9) and (I.7.10). From (I.7.13) and (I.7.9) squared one may write E' = (1/2)I'12 + { (1/2)I'1+ Mgrcmscosθ } or α = 2 + {+ βcosθ } = 2 + Ve(θ) α = 2E'/I'1, β = 2Mgrcms/I'1 (I.7.18) where Ve(θ) is the effective potential for the top problem. We can then carry out the same qualitative turning point analysis as shown in Fig (C.5.17) and we find that θ(t) bounces back and forth between two angles θ1 and θ2. The tip of the top then has these typical motion patterns depending on the size of the four parameters a,b,α,β (GPS p 215) : (I.7.19) The pattern on the left is usually called nutation during precession. In the middle pattern, changes sign during the run from θ1 to θ2. See Thornton and Marion pp 456-460 for details of fast and slow precession and of the above patterns. Reader Exercise: Solve this problem directly from the equations of motion (I.4.3) where (N)'3 = 0. Here is a comparison of characteristics of the spinning top and the torque-free rotor of Section I.5 ; Torque-free rotor Spinning Gravity Top (ω)'3 = constant = ωcosα (ω)'3 = constant = (aI'1/I'3) (ω)'1 = ωsinα cos(Ω't) (ω)'1 = sinθsinψ + cosψ (ω)'2 = ωsinα sin(Ω't) (ω)'2 = sinθcosψ - sinψ ω2 = constant ω2(t) = 2 sin2θ + 2 + (aI'1/I'3)2 θ = constant = tan-1[ (I'1/I'3) tanα ] θ(t) = involves elliptic functions = constant = - Ω' (t) = a (I'1/I'3) - cosθ( b - cosθ a)/sin2θ = constant = (L/I'1) (t) = T = constant T + Mgrcmscosθ = constant L = constant Lz = constant = aI'1 Lz' = constant = bI'1 ω1(t) = ω sin(θ-α) cos[(L/I'1)t -π/2] ω1(t) = sinθsinφ + cosφ ω2(t) = ω sin(θ-α) sin[(L/I'1)t -π/2] ω2(t) = - sinθcosφ + sinφ ω3 = ω cos(θ-α) ω3(t) = cosθ + . (I.7.20) I.8 Gravitational Torque on an oblate Earth: Precession of the Equinoxes For motivation, we consider the gravitational torque due to a large distant mass M on an circle of equally-spaced point masses m (only 4 are shown), with everything lying in a plane The masses m are glued to the rigid circle and we ignore their gravitational attraction to each other. The large distance from mass M to the circle center is b. (I.8.1) We consider the mass pair AB to be a "dumbbell" of masses, and we have already computed in (F.3.13) the torque (about the circle center) on this dumbbell due to the gravity of mass M, assuming b >> r, NAB = - 6GMmb-3r2sinθcosθ . // μ2 = 1/2 (F.3.13) (I.8.2) Because the mass A is closer to M than the mass B, the torque on mass A going into the plane of paper (right hand rule N = r x F ) is larger than the torque on mass B going out of the plane of paper, and the resultant torque on the dumbbell AB is into the plane of paper. This torque wants to restore the dumbbell AB to an aligned orientation, just as it does for any dumbbell satellite. However, the torque on the dumbbell CD is just the opposite of the torque on AB : NAB = - 6GMmb-3r2sinθcosθ NCD = - 6GMmb-3r2sin[-θ]cos[-θ] = = + 6GMmb-3r2sinθcosθ so NAB + NCD = 0 . (I.8.3) When we sum over all pairs of masses on the circle, the net result is that mass M exerts zero torque on the circle of masses. Consider now a tilted ellipse of equal point masses: (I.8.4) We now have NAB = - 6GMmb-3r12sinθcosθ NCD = - 6GMmb-3r22sin[-θ]cos[-θ] = = + 6GMmb-3r22sinθcosθ NAB + NCD = 6GMmb-3(r22 - r12) sinθcosθ . (I.8.5) Since r2 > r1, the total torque exerted by M on the pairs AB and CD is non-zero and is directed toward the viewer. When we sum over all pairs of masses on the ellipse, the net result is that mass M exerts a torque toward the viewer which wants to restore the ellipse of masses to an aligned position with the long ellipse axis lying along the x axis. With this planar example as motivation, consider the following 3D picture of gravitational mass M orbiting the oblate red Earth on a blue orbit (circular or elliptical), (I.8.6) At the instant shown in the drawing, mass M (imagined to lie on the Celestial Sphere) is at some latitude δ (declination) relative to the Earth's equatorial plane. Azimuthally, M is located at longitude α (right ascension) relative to some specified meridian (vernal equinox) on the Celestial Sphere. The unit vector is the usual azimuthal unit vector which appears as in Fig (E.1.1). In particular, = (-sinα,cosα,0) as shown in (E.2.6). Based on our 2D example, we expect mass M to exert a torque on the Earth in the sense of the black arrows. This torque is attempting to rotate the Earth toward a position where M then lies on the Earth's equatorial plane. If M is the Sun then, as the Sun orbits, the declination δ moves between +23.44o and -23.44o in our current era. If M is the Moon, then, as the Moon orbits, declination δ moves between 28.58o and -28.58o near the major lunar standstill, and between 18.30o and -18.30o near the minor lunar standstill which occurs 9.3 years after the major one. See Fig (8.8.58a) and (8.8.58b) which show these two standstills (the wiki lunar standstill page is very good on this subject). The range variation is caused by the fact that the Moon's orbital plane is tilted 5.14o relative to the Sun's orbital plane causing the 23.44 ± 5.14 ranges. The torque N of mass M on the Earth is stated in Williams equation (2) to be, N = - 3GM b-3(I3 - I1) sinδ cosδ = . (I.8.7) The parameters here correspond go Fig (I.8.6) and I3 > I1 are the two moments of inertia for the oblate Earth (which of course is axisymmetric). The result is valid for any axially symmetric mass distribution (including prolate). If δ is replaced by polar angle θ = π/2 - δ, one can replace, sinδ cosδ = sin(π/2-θ) cosδ(π/2-θ) = cosθ sinθ = sinθ cosθ . We shall derive (I.8.7) in the next section. For now, one can see that for α = π/2 one will have = - so N = (positive) which is a torque trying to rotate the Earth CCW to become aligned with the mass M line (length b). In this case Nx = 3GM b-3(I3 - I1) sinθcosθ which bears a resemblance to our toy planar example result (I.8.5) that (NAB + NCD) = 6GMb-3(mr22 - mr12) sinθcosθ . Precession of the Equinoxes and the Pope We think now of Frame S as an inertial frame associated with the Sun, and rotating Frame S' being embedded in the Earth. When William's torque is put into the equations of motion (I.4.3) with I'1 = I'2, one finds that the Earth's ω vector would precess relative to the stars with a period of about 81,000 years. But the Moon also orbits the oblate Earth, and causes the Earth to precess as well. As with the tides, the Moon's effect on precession is about twice that of the Sun's. The net result of the Sun and the Moon is that the Earth's rotation axis precesses with a period of around 25,800 years. These precession periods are long because the Earth is only slightly oblate and the resulting gravitational torques are relatively small. The net effect is known as "axial precession" and historically as "the precession of the equinoxes" since the seasons of the year slide about one degree per 72 years (~ 25,800/365), or one day every 129.4 years. Since this slippage threatened to pull Easter away from the vernal equinox (Spring), Pope Gregory in 1582 set up the Gregorian calendar to neutralize the slippage. In this scheme, one has no leap year every XX00 year except when XX is a multiple of 4. This knocks out 3 days every 400 years, and 400/3 = 133.3 ≈ 129.4. I.9 Derivation of the Oblate Earth torque formula In this Section we take a whirlwind tour of "advanced" potential theory and then abstract out the relatively simple pieces needed for our problem. Sturm Liouville Transforms Away from sources of charge, the electrostatic potential must satisfy the Laplace equation. The same is true for the gravitational potential away from mass sources. In general any such potential can be expressed as a sum over the "atomic forms" (harmonics) of a given coordinate system. For spherical coordinates the Laplace equation atomic forms can be taken as [rn , r-n-1] [ Pnm(z), Qnm(z) ] [eimφ, e-imφ] z = cosθ (I.9.1) where n = 0,1.2.. and m = -n,-n+1,..0,1,..n. These forms are appropriate for a region of space which has the full range of angle φ. In such a region, for the potential to be single valued in azimuth one must have m = integer so eim(φ+2π) = eimφ. For technical reasons, this in turn forces n to be an integer n ≤ m. If θ = 0 is part of the range of interest, one cannot have Qnm(z) terms since Qnm(z) blows up at z = 1. In this case one has a reduced set of atoms to work with, [rn , r-n-1] [ Pnm(z)] [eimφ, e-imφ] z = cosθ . (I.9.2) The associated Legendre functions Pnm(z) are solutions to a standard-issue Sturm Liouville problem and as such one can write down a transform (expansion and projection), an orthogonality condition, and a completeness condition. g(z) = Σn=|m|∞gnm Pnm(z) z = cosθ // expansion gnm = (1/Knm)!Syntax Error, Idz Pnm(z) g(z) Knm = (n+1/2)-1 f(n,m) // projection !Syntax Error, Idz Pnm(z)Pn'm(z) = δn,n' Knm n,n' = |m|, |m|+1, ... // orthogonality Σn=|m|∞(1/Knm) Pnm(z') Pnm(z) = δ(z'-z) f(n,m) = Γ(n+m+1)/Γ(n-m+1) . // completeness (I.9.3) For verification of orthogonality, see Jackson (3.52). For example, one may expand g(z) onto the Pnm(z) with coefficients gnm which in turn are given by the projection on the second line. The various constants are shown on the right. In these equations, the parameter m is somewhat of a bystander parameter. For m = 0 one sees that f(n,0) = 1 and then Kn0 = 1/(n+1/2) so setting gn0 = gn one gets, g(z) = Σn=0∞gn Pn(z) // expansion gn = (n+1/2)!Syntax Error, Idz Pn(z) g(z) // projection !Syntax Error, Idz Pn(z)Pn'(z) = δn,n'/ (n+1/2) // orthogonality Σn=0∞(n+1/2) Pn(z') Pn(z) = δ(z'-z) . // completeness (I.9.4) Expression for 1/R in Spherical Coordinates In potential theory, if one places a point source q at location r' and views it from point r, the potential at r is given by V(r) = k R = |r-r'| . (I.9.5) For this potential, we define two regions of space: inside and outside of the " r' circle" , as shown on the left below. (I.9.6) In each of these regions 1/R, being a solution of the Laplace equation, must be expressible as a linear combination of the atomic forms (the sum is Σmn = Σn=0∞ Σm=-nn) , 1/R = Σmn Amn(r',θ',φ') rn Pnm(cosθ) eimφ for r inside the r' circle 1/R = Σmn Bmn(r',θ',φ') r-n-1 Pnm(cosθ) eimφ for r outside the r' circle . (I.9.7) Inside the r' circle the powers r-n–1 for n = 0,1,2.. are ruled out since they blow up at r = 0. Outside the r' circle the powers rn for n = 1,2.. are ruled out since they blow up as r→ ∞. Since R is completely symmetric under r ↔ r' one must also be able to write, 1/R = Σmn Amn(r,θ,φ) r'n Pnm(cosθ') e-imφ' for r' inside the r circle 1/R = Σmn Bmn(r,θ,φ) r'-n-1 Pnm(cosθ') e-imφ' for r' outside the r circle . (I.9.8) Here we have taken the liberty to change eimφ' to e-imφ'since we know that the imaginary parts make no contribution to 1/R since R is real and functions and coefficients are real. The implication is that one must be able to write 1/R = Σmn Cmn (rn Pnm(cosθ) eimφ)( r'-n-1 Pnm(cosθ') e-imφ') for r inside the r' circle 1/R = Σmn Dmn ( r-n-1 Pnm(cosθ) eimφ)(r'n Pnm(cosθ') e-imφ') for r outside the r' circle or (I.9.9) 1/R = Σmn Cmn rn r'-n-1 Pnm(cosθ)Pnm(cosθ')eim(φ-φ') for r inside the r' circle 1/R = Σmn Dmn r'n r-n-1Pnm(cosθ)Pnm(cosθ')eim(φ-φ') for r outside the r' circle . As r'→r, one sees that the two coefficient sets must be the same, Cmn= Dmn. The final result is then usually written as a single line in this obvious manner, where r> = max(r,r') and r< = min(r,r'), 1/R = Σmn Cmn (r<)n (r>)'-n-1 Pnm(cosθ)Pnm(cosθ')eim(φ-φ') . (I.9.10) By using a small and thin Gaussian box in the right place, throwing around some delta functions and curvilinear scale factors, and making appropriate invocations about jump conditions across the box, one can evaluate the coefficients Cmn in the equivalent of the above 1/R atomic expansion for any curvilinear coordinate system. For spherical coordinates one finds, Cmn = f(n,-m) ≡ Γ(n-m+1)/Γ(n+m+1) . (I.9.11) Thus 1/R = Σn=0∞ Σm=-nn f(n,-m) (r<)n (r>)'-n-1 Pnm(cosθ)Pnm(cosθ')eim(φ-φ') . (I.9.12) By reflecting the negative m values to the positive side, this may also be expressed as 1/R = Σn=0∞ Σm=0n εm f(n,-m) (r<)n(r>)-n-1 Pnm(z) Pnm(z') cos[m(φ-φ')] (I.9.13) where εm = 2 - δm,0 is the so-called Neumann factor (ε0 = 1, εother = 2) . Potential of a Source Distribution In potential theory, if there is some isolated source distribution ρ(r) (charge in electrostatics, mass in gravitation), the potential can be written as an integral over the source distribution in this manner, V(r) = k ∫dV' = k ∫dV' . (I.9.14) There are no boundary conditions, the source sits in empty infinite space. This is just the superposition of V = k(q/R) for a point source q where q = ρ(r')dV' . In electrostatics the constant k is 1 in cgs units, and is 1/4πε0 in SI units. For gravity k = -G. Inserting the 1/R expansion above we find that for r outside the r' circle. V(r) = k ∫dV' ρ(r') {Σn=0∞ Σm=-nn f(n,-m) (r<)n(r>)-n-1 Pnm(cosθ)Pnm(cosθ')eim(φ-φ')} = k Σn=0∞ Σm=-nn f(n,-m) Pnm(cosθ)eimφ r-n-1∫dV' ρ(r') r' n Pnm(cosθ') e-imφ' . (I.9.15) We note in passing that one usually defines the "spherical harmonics" in this manner Ynm(θ,φ) ≡ Pnm(cosθ)eimφ // Jackson (3.53) (I.9.16) Therefore, f(n,-m) Pnm(cosθ)eimφ Pnm(cosθ') e-imφ' = Ynm(θ,φ) Y*nm(θ',φ') (I.9.17) and then V(r) = k Σn=0∞ Σm=-nn f(n,-m) Ynm(θ,φ) r-n-1∫dV' ρ(r') r' n Y*nm(θ',φ') (I.9.18) 1/R = Σn=0∞ Σm=-nn (r<)n (r>)'-n-1 Ynm(θ,φ) Y*nm(θ',φ') . // Jackson (3.70) (I.9.19) One can define the exterior "multipole moments" Tnm as Tnm ≡ ∫dV' ρ(r') r' n Y*nm(θ',φ') // dV' = r'2dr' sinθ'dθ'dφ' (I.9.20) and then the V(r) expansion becomes, V(r) = k Σn=0∞ Σm=-nn f(n,-m) Ynm(θ,φ) r-n-1 Tnm (I.9.21) The interior multipole moment expansion then has rn ↔ r-n-1. Although for a given n there are 2n+1 multipole moments, the historical names for n = 0,1,2,3 are monopole, dipole, quadrupole and octupole. These terms refer to simple source patterns in Cartesian space which have such moments. We will be dealing with n = 2 quadrupole below. For a mass distribution (the oblate Earth) which is azimuthally symmetric, only the Tn0 moments are non-vanishing due to the e-imφdφ integration in (I.9.20). The exterior expansion (I.9.15) then becomes, V(r) = k Σn=0∞ Pn(cosθ) r-n-1∫dV' ρ(r',θ') r' n Pn(cosθ') // f(n,0) = 1 (I.9.22) and we define some new moments Jn such that (Jn are constants, not Bessel functions) Jn ≡ k∫dV' ρ(r') r' n Pn(cosθ') // projection V(r) = Σn=0,1.. Pn(cosθ) r-n-1 Jn . // expansion (I.9.23) If the source distribution ρ(r') is invariant under vertical reflection, as is the case for the oblate Earth, the odd Jn vanish. For example, this knocks out the n = 1 dipole term. The decay factor r-n-1 becomes more severe for larger n, so "far away" one takes only the first few terms in the expansion. We shall take the first two contributing terms, V(r) ≈ P0(cosθ) r-1 J0 + P2(cosθ) r-3 J2 . // P0(z) = 1 (I.9.24) The moments here are J0 = -G ∫dV' ρ(r') r' 0 P0(cosθ') = -G ∫dV' ρ(r',θ') = -GME J2 = -G∫dV' ρ(r') r' 2 P2(cosθ') (I.9.25) so V(r) ≈ -GME/r + J2 P2(cosθ)/r3  . (I.9.26) The first term is the potential of the spherical Earth, while the second term is a correction. Our next task is to compute the moment J2. Calculation of the moment J2 We shall assume for our azimuthal mass distribution that any line segment from the origin to a point on the boundary lies within the boundary. The region is "star-like" since all star rays from the origin lie in the region. In this case, when we do the volume integration to compute moments, the θ and φ integrations have full range, but the r integration runs from r = 0 to some r = r(θ) which defines the boundary. This boundary allows for both oblate or prolate shape types. Removing the dummy primes in (I.9.25) one has, J2 = -G ∫dV ρ(r) r2 P2(cosθ) = -2πG !Syntax Error, Idθ sinθ P2(cosθ)!Syntax Error, Idr r4 ρ(r,θ) . (I.9.27) We now digress to recall some facts about moments of inertia from Section I.1 above. In particular, Iij ≡ ∫dx1dx2dx3 ρ(x)[ r2δij - xixj] = ∫dV ρ(r)( r2δij - xixj) . (I.1.3) In particular, since our "planet" is axisymmetric, the inertia tensor is diagonal with I1 = I2 and I33 = I3 = ∫dV ρ(r)( r2 - z2) I11 = I1 = ∫dV ρ(r)( r2 - x2) = I22 = I2 . (I.9.28) Therefore, I3 - I1 = ∫dV ρ(r)(x2-z2) = ∫dV ρ(r)([rsinθcosφ]2-[rcosθ]2) = ∫dV ρ(r) r2(sin2θ cos2φ - cos2θ) = !Syntax Error, Idφ !Syntax Error, Idθ sinθ (sin2θ cos2φ - cos2θ)!Syntax Error, Ir4dr ρ(r,θ) = !Syntax Error, Idθ sinθ (π sin2θ - 2πcos2θ)!Syntax Error, Ir4dr ρ(r,θ) // !Syntax Error, Idφ cos2φ = π = π!Syntax Error, Idθ sinθ (sin2θ - 2cos2θ)!Syntax Error, Ir4dr ρ(r,θ) = π!Syntax Error, Idθ sinθ (1 - 3cos2θ)!Syntax Error, Ir4dr ρ(r,θ) = -2π!Syntax Error, Idθ sinθ P2(cosθ)!Syntax Error, Ir4dr ρ(r,θ) . (I.9.29) Comparing this with the J2 integral in (I.9.27) one sees that for any azimuthally symmetric mass distribution, J2 = G(I3-I1) (I.9.30) and then the potential is (if also vertically mirror symmetric so no P1 term), V(r) ≈ - GME/r + G (I3-I1) P2(cosθ)/r3. (I.9.31) This result is valid for any "star" perimeter r(θ) and for any mass distribution ρ(r,θ). We shall come back later and compute I3-I1 for a slightly oblate or prolate spheroid of uniform density. But first we want to derive the torque equation (I.8.7) as promised. Calculation of the Torque In Fig (I.8.6) we show mass M which is "far away" from the Earth. The Earth generates a potential V(r) as computed in (I.9.31). This potential creates a force F = -MV(r) on mass M. That force in turn creates a torque r x F acting on mass M. Mass M of course then creates a reverse torque N = - r x F on the Earth. We first compute V(r)as follows: V(r) = -GME(1/r) + G (I3-I1) [ ( P2(cosθ) )/r3 + P2(cosθ)(1/r3) ] Since f(r) = f'(r), the terms with (1/r) and (1/r3) have no effect on the torque. One then has, N = - r x F = M r x V(r) = GM (I3-I1) r-3 r x (P2(cosθ)) . (I.9.32) In spherical coordinates, ( P2(cosθ) ) = (1/r) ∂θP2(cosθ) . (I.9.33) Then, N = - r x F = GM (I3-I1) r-3 r x [ (1/r) ∂θP2(cosθ) = GM (I3-I1)r-3 ∂θP2(cosθ) x = GM (I3-I1)r-3 ∂θP2(cosθ) // (E.2.12) = GM (I3-I1)r-3 ∂θ{(1/2)(3cos2θ - 1)} = GM (I3-I1)r-3 (3/2)∂θ{cos2θ} = GM (I3-I1)r-3 (3/2)[-2cosθ sinθ ] = - 3GM r-3 (I3-I1) cosθ sinθ . (I.9.34) This is the result we quoted from Williams in (I.8.7) with r = b and = . Description of the perimeter of a slightly oblate or prolate spheroid We now digress a bit on ellipses and spheroids created from them. Consider an ellipse of the form x2/A2 + z2/C2 = 1, (I.9.35) If we measure polar angle θ down from the z axis, then x = rsinθ and z = rcosθ and the ellipse equation becomes r2 = . (I.9.36) For an "oblate" (wide) ellipse we have A > C. If the ellipse is only slightly oblate, write A = C + δ . For oblate δ > 0, and for prolate δ < 0. In this case one gets, r2 = ≈ ≈ ≈ A2 [ 1-2(δ/C)cos2θ ] ≈ A2 [ 1-2(δ/A)cos2θ ] . Then r ≈ A [ 1-(δ/A)cos2θ ] . (I.9.37) The "eccentricity" e of the ellipse is e2 ≡ (A2-C2)/A2 = (A+C)(A-C)/A2 ≈ 2Aδ/A2 = 2δ/A . One then has r = A [ 1- (1/2)(2δ/A)cos2θ ] = [ 1- (1/2)e2cos2θ ] . On the other hand, the "ellipticity" of the ellipse is (also known as the flattening parameter) ε ≡ (A-C)/A = δ/A // ε = e2/2 (I.9.38) where ε > 0 for oblate and ε < 0 for prolate. Then, r(θ) = A [ 1- (1/2)(2ε)cos2θ ] = A[ 1- ε cos2θ ] . (I.9.39) We now form a slightly oblate spheroid by rotating our slightly oblate ellipse about the vertical axis, introducing an azimuthal coordinate φ. The average value of r over this spheroid is given by, a ≡ <r> = A[ 1- ε <cos2θ> ] . (I.9.40) At this point assume a uniform mass density ρ(r) = ρ so that <cos2θ> = = = ≈ = = = 1/3 (I.9.41) where corrections of order ε can be ignored since these create order ε2 in a. Then from (I.9.40), a ≡ <r> = A[ 1- ε (1/3) ] . (I.9.42) Now process r(θ) of (I.9.39) as follows, r(θ) = A[ 1- ε cos2θ ] = [ 1- ε cos2θ ] ≈ a(1+ε/3)(1- ε cos2θ) ≈ a [ 1 + ε/3 - ε cos2θ ] = a [ 1 - (ε/3)(3cos2θ - 1) ] = a [ 1 - (2/3)ε P2(cosθ) ] // since P2(cosθ) = (1/2)(3cos2θ - 1) or r(θ) = a [ 1 - (2/3)ε P2(cosθ) ] . (I.9.43) One can regard r(θ) as a function defining the perimeter of the slightly oblate or prolate spheroid of ellipticity (flatness) ε. If ε = 0, one of course finds that r = a for a sphere. The reader seeking verification will find (I.9.43) to be the very first equation on this nice web page http://farside.ph.utexas.edu/teaching/336k/Newtonhtml/node108.html which we have used as a guide throughout. These are notes by Richard Fitzpatrick. Calculation of I3- I1 Recall from (I.9.29) that I3 - I1 = -2π!Syntax Error, Idθ sinθ P2(cosθ)!Syntax Error, Ir4dr ρ(r,θ) = -2πρ!Syntax Error, Idθ sinθ P2(cosθ) (1/5) r(θ)5 . (I.9.44) Now r(θ) = a [ 1 - (2/3)ε P2(cosθ) ] so r(θ)5 ≡ a5 [ 1 - 5 (2/3)ε P2(cosθ) ] . (I.9.45) Then I3 - I1 = -2π ρ (1/5) a5!Syntax Error, Idθ sinθ P2(cosθ) [ P0(cosθ) - 5 (2/3)ε P2(cosθ) ] . (I.9.46) Recall from (I.9.4) that !Syntax Error, Idz Pn(z)Pn'(z) = δn,n'/ (n+1/2) . (I.9.4) The first integral in (I.9.46) thus vanishes leaving I3 - I1 = -2πρ (1/5) a5 [ - 5 (2/3)ε] {!Syntax Error, Idθ sinθ P2(cosθ) P2(cosθ) } = 2πρ (1/5) a5 5 (2/3)ε {1/[5/2]} = 2πρ (1/5) a5 5 (2/3)ε (2/5) = 2πρ (1/5) a5 (2/3)ε (2) = (8π/15)ερa5 . (I.9.47) To zeroth order the mass density is that of a sphere, so ρ = ME/V = ME / [ (4/3)πa3] (I.9.48) and then I3 - I1 = (8π/15)εa5 * ME / [ (4/3)πa3] = (8/15)εa2 * ME / [ (4/3)] = (8/15)εa2 * ME (3/4) = (2/5) ε MEa2 (I.9.49) and the moment J2 is then J2 = G(I3 - I1) = (2/5) ε GMEa2 . (I.9.50) The potential can then be written V(r) = -GME/r + G (I3-I1) P2(cosθ)/r3. = -GME/r + (2/5) ε GMEa2 P2(cosθ)/r3 (I.9.51) Reader Exercise: Show that a uniform sphere of radius a and mass M has a moment of inertia about any axis of I = (2/5) Ma2 . For the Earth then one has roughly I1 ≈ I3 ≈ (2/5) Ma2 and then (I.9.49) states, I3 - I1 = (2/5) ε MEa2 ≈ ε I1 so ε ≈ // "McCullough's formula" (I.9.52) which relates the ellipticity to the moments of inertia for an axially symmetric "planet". Reader Exercises: 1. Compute the integrals I1 and I3 separately as in (I.9.29) and conclude that McCullough's formula is valid for any "star-like" shape and any azimuthally symmetric mass density distribution ρ(r,θ). 2. Show that McCullough's formula is valid even if the boundary is not "star-like". 3. For a uniform arbitrary oblate spheroid (not one that is just slightly oblate) compute, with no approximations, the exact external gravitational potential and resulting torque exerted by an external point mass M. Use oblate spheroidal coordinates. A candidate potential solution appears on page 62 of MacMillan (which may be available online). That solution is quoted here, http://scienceworld.wolfram.com/physics/OblateSpheroidGravitationalPotential.html In the electrostatics problem of a charged conducting spheroid, the mobile source (charge) is distributed on the spheroidal surface, while in the gravitational problem the fixed mass is distributed throughout the spheroidal volume. How does this fact affect the two problem solutions?