Historical and Modern Perspectives on Hamilton-Jacobi Equations

Historical and Modern Perspectives on Hamilton-Jacobi Equations

Robin Ekman

June 10, 2012

Abstract

We present the Hamilton-Jacobi equation as originally derived by Hamilton in 1834 and 1835 and its modern interpretation as determining a canonical trans- formation. We show that the method of characteristics for partial differential equations, applied to the Hamilton-Jacobi equation, yields Hamilton’s canonical equations and the action as the solution. Canonical perturbation theory is ap- plied to the Sun-Earth-Jupiter system to calculate perturbations of the orbital elements. We also study tidal locking and find a damped harmonic oscillator equation for a moon that is almost locked. Finally, we present a modern (1980s), weaker notion of solutions to the Hamilton-Jacobi equation, viscosity solutions.

We reproduce proofs that this concept is consistent with the classical concept

and that visocsity solutions are unique, for locally Lipschitzian Hamiltonians.

Contents

I Historical Background for Hamilton-Jacobi Equations in Clas-

sical Mechanics 4

I.1 Hamilton’s Characteristic and Principal Functions . . . . 4

I.2 The Principal Function as a Canonical Transformation . . . . 10

I.3 Hamilton’s Equations as Characteristics . . . . 13

II Classical Perturbation Theory Applied to Planetary Orbits 15 II.1 Canonical Perturbation Theory . . . . 15

II.2 Separability and Action-Angle Variables . . . . 17

II.3 Action-Angle Variables for the Kepler Problem . . . . 18

II.4 Perturbation . . . . 24

II.5 Averaging . . . . 25

II.6 Perturbations of Orbital Elements . . . . 26

II.7 Tidal Locking . . . . 29

III Existence and Uniqueness of Solutions 32 III.1 Viscosity solutions . . . . 33

III.2 Classical Solutions are Viscosity Solutions . . . . 34

III.3 Viscosity Solutions are Unique . . . . 36

III.4 Other Methods . . . . 38

A Expansion of the Perturbation of the Inclination 39 List of Figures II.1 Orbital elements . . . . 22

II.2 Angles used for the Kepler problem. ON is the line of nodes; R is the position of the orbiting body. . . . 22

II.3 Eccentric anomaly ψ and true anomaly ϕ. . . . . 23

II.4 Average precession of the perihelion using values for the Sun,

Earth and Jupiter. . . . 27

II.5 Trajectory for a planet perturbed by another planet . . . . 27

II.6 An orbiting body of finite extent. (Not to scale.) . . . . 30

III.1 A solution to the wave equation? . . . . 32

[Lagrange showed] that the most varied consequences respecting the motions of systems of bodies may be derived from one radical for- mula; the beauty of the method so suiting the dignity of the results, as to make of his great work a kind of scientific poem. But the science of force, or of power acting by law in space and time, has undergone already another revolution, and has become already more dynamic ... and while the science is advancing thus in one direction by the improvement of physical views, it may advance in another direction also by the invention of mathematical methods.

—William Rowan Hamilton, On a General Method in Dynamics, 1834.

Chapter I

Historical Background for Hamilton-Jacobi Equations in Classical Mechanics

In this chapter we will consider three perspectives on Hamilton-Jacobi equations in classical mechanics. The common denominator in these three presentations is showing that a system of ordinary differential equations (ODE) can be reduced to one partial differential equation (PDE), and vice versa. We shall discuss mechanical problems by which we mean a system of particles, interacting in some way described by a potential function. By a solution we mean a set of functions of time giving the motion of each particle, given the initial positions and velocities.

First we will present Hamilton’s original derivation from his 1834 and 1835 papers On a General Method in Dynamics and A Second Essay on a General Method in Dynamics [9, 10] using the calculus of variations, extending the work of Lagrange. An important step here is Hamilton’s derivation of a system of first-order ODE, the canonical equations. Secondly, we present a later view where the Hamilton-Jacobi equation determines a coordinate transformation that also solves the mechanical problem at hand [7, 1]. This view assumes the canonical equations to have been found but uses no more methods from the calculus of variations. Lastly, we show that if we are given the Hamilton- Jacobi equation, the method of characteristics in the theory of PDE generates Hamilton’s canonical equations [6].

I.1 Hamilton’s Characteristic and Principal Func- tions

Hamilton’s General Method

In his 1834 paper [9] Hamilton considers a system of n particles (points, in his

terminology) interacting by repulsion or attraction depending on their distances.

That is, from the equations of motion

n

X

i=1

m i x ⁰⁰ _i · δx i = δU

where m i denote the masses, x ⁰⁰ _i is the acceleration in rectangular (Cartesian) coordinates (likewise, x ⁰ _i will be the velocity), δx i an infinitesimal displacement and δU an infinitesimal change in the function U , which describes the interac- tions as

U = X

m i m j f (r ij ).

We have, unlike Hamilton, used vector notation and the dot product, and will continue to do so, as the modern reader should be familiar with this condensing notation. Hamilton called U the ‘force-function’, but we know it today as the potential and conventionally define it with the opposite sign; we would write mx ⁰⁰ = −∇U. We shall follow Hamilton’s choice of sign.

Next Hamilton considers

T = 1 2

X m i x ⁰ _i ² ,

x ⁰ _i being the velocity, which he calls half “the living force of the system”, in modern terms the kinetic energy. Then using the equations of motions we can write

dT = dT

dt dt = X

m i x ⁰ _i · x ⁰⁰ i dt = X

m i dx i · x ⁰⁰ i = dU and we can conclude that

T = U + H (I.1)

where H is a constant of the motion. Hamilton calls this the “law of living force”; remembering that his definition of U differs from the modern by a sign we recognize it as the conservation of energy.

If the initial conditions change, then H can change too and we have δT = δU + δH. Then integrating along the path taken by the system

Z X

m i dx i · δx ⁰ i =

Z X

m i dx ⁰ _i · δx i + Z

δHdt. (I.2)

Hamilton presents his next result as following from “the principles of the calculus of variations”, skipping some steps. If we define, as Hamilton does,

V = Z t

0

2T dt = Z

X m i x ⁰ _i · dx i (I.3)

the variation is

δV =

Z X

m i δx ⁰ _i · dx i + X

m i x ⁰ _i · δdx i .

The second term can be integrated by parts, putting another derivative on x:

δV =

Z X

m i δx ⁰ _i · dx − Z t

0

X m i x ⁰⁰ _i · δx i dt +  X m i x ⁰ _i · δx i  x

a

where we also used that the order of variation and differentiation can be inter- changed. Now the first integral can be replaced using (I.2) and since H is a constant of motion we can replace the integral with H with tδH. That is,

δV =

Z X

m i δx i · dx ⁰ i + tδH − Z t

0

X m i x ⁰⁰ _i · δx ⁱ dt +  X m i x ⁰ _i · δx ⁱ  x

a

. But the two integrals cancel, since x ⁰⁰ _i dt = dx ⁰ _i and we are left with

δV = X

m i x ⁰ _i · δx i − X

m i a ⁰ _i · δa i + tδH

where a i is the initial position. If V is considered a function of the initial and end points, we have the system of equations

∂V

∂x i

= m i x ⁰ _i (I.4)

∂V

∂a i = −m ⁱ a ⁰ _i (I.5)

∂V

∂H = t (I.6)

which describes the system in terms of one single function, the characteristic function V .

Now using this system of equations in the law of living force (I.1) Hamilton obtains two PDE, for V , the former which he calls final and the latter initial

1 2

X 1

m i

 ∂V

∂x i

 2

= U + H (I.7)

1 2

X 1

m i

 ∂V

∂a i

 2

= U 0 + H (I.8)

where U 0 is the inital force function. If the function V were known, the me- chanical problem would be solved, for we can verify that the system of equations for the characteristic function generates the same equations of motions that we started with. Then the problem of solving n ODE of second order has been transformed to solving one first order PDE.

In the last section of his first paper, Hamilton introduces the function S defined by the relation

V = tH + S.

Then S, which Hamilton here calls the auxiliary function, using (I.1), (I.3) and tH = R t

0 H dt since H is a constant of motion, satisfies S =

Z t 0

(T + U )dt (I.9)

and its partial derivative with respect to t is ^∂S _∂t = −H, its other partial deriva-

tives satisfy the same equations as do the partial derivatives of V .

Hamilton’s Second Essay

In the Second Essay, Hamilton begins by showing that if the 3n Cartesian coordinates are “functions of 3n other and more general marks of position”, η 1 , η 2 , . . . , η 3n , the form of the equations of motion is that given by Lagrange

d dt

∂T

∂η _i ⁰ − ∂T

∂η i

= ∂U

∂η i

(I.10) where T is now a function of the η i and η _i ⁰ . This is a set of n ODE of second order.

Since for each x ⁰ _i we have x ⁰ _i = P η _j ^{0 ∂x} _∂η

j

(Hamilton did not consider time- dependent coordinate systems) and T is homogeneously quadratic when ex- pressed in x ⁰ _i , it is homogeneously quadratic in η _j ⁰ and a theorem of Euler ¹ gives

2T = X η ⁰ _j ∂T

∂η ⁰ _j . (I.11)

Hamilton then notes that taking the variation of this, and subtracting the vari- ation of T in the form

δT = X  ∂T

∂η _i ⁰ δη ⁰ _i + ∂T

∂η i δη i



(I.12) one obtains

δT = X  η ⁰ _i δ ∂T

∂η ⁰ _i − ∂T

∂η i

δη i



. (I.13)

In the following paragraph Hamilton writes

“If then we put, for abridgement,

∂T

∂η ₁ ⁰ = ω 1 , . . . ∂T

∂η ⁰ _3n = ω 3n

and consider T (as we may) as a function of the following form, T = F (ω 1 , ω 2 , . . . ω 3n , η 1 , η 2 , . . . η 3n ),

”

and the canonical formalism is introduced. Using these variables the relations

∂F

∂ω i = η _i ⁰ ∂F

∂η i = − ∂T

∂η i

follow from taking the variation of F : δF = X ∂F

∂ω i

δω i + ∂F

∂η i

δη i = X

η _i ⁰ δω i − ∂T

∂η i

δη i

and comparing with the variation of T , (I.13). Lagrange’s equations of motion (I.10) become

dω i

dt = ∂(U − F )

∂η i

. (I.14)

1

The theorem states: if f (x) : R

→ R is homogeneous of degree n (f (tx) = t

f (x) for all

t ∈ R) and has continuous derivatives, then nf (x) = x ·

.

Now putting H = F − U we have two sets of first order ODE dη i

dt = ∂H

∂ω i

dω i

dt = − ∂H

∂η i

.

These are now known as Hamilton’s (canonical) equations, and we usually write q instead of η and p instead of ω.

Hamilton now defines the principal function S by the integral S =

Z t 0

X ω i

∂H

∂ω i − H dt (I.15)

which using the canonical equation can be seen to be equal to his previous definition (I.9). Hamilton takes the variation of S, keeping t and dt constant so that

δS = Z t

0

δS ⁰ dt where S ⁰ is the integrand, thus

δS ⁰ = X 

ω i (δ ∂H

∂ω i ) − ∂H

∂η i

δη i + ∂H

∂ω i

δω i − ∂H

∂ω i

δω i

 , which using Hamilton’s equations is

δS ⁰ = X  ω i δ dη i

dt + dω i

dt δη i



= d dt

X ω i δη i

so that the variation of its integral is δS = X

ω i δη i − p i δe i

where p i and e i are initial values of ω i and η i , respectively. Hamilton thus concludes that the principal function S satisfies

ω i = ∂S

∂η i

p i = − ∂S

∂e i

. (I.16)

Next, Hamilton computes ^∂S _∂t from dS

dt = ∂S

∂t + X ∂S

∂η i

dη i

dt

which using Hamilton’s canonical equations and the definition of S gives

∂S

∂t = S ⁰ − X ω i

∂H

∂ω i = −H.

Hamilton here notes that −H is a constant of motion, as he considers only conservative systems. It was noted by Jacobi [11] that this restriction is not necessary.

Hamilton concludes that the principal function S must satisfy a pair of PDE

∂S

∂t + F  ∂S

∂η 1

, ∂S

∂η 2

, . . . ∂S

∂η 3n

, η 1 , η 2 , . . . η 3n



= U (η 1 , η 2 , . . . η 3n ) (I.17)

∂S

∂t + F  ∂S

∂e 1

, ∂S

∂e 2

, . . . ∂S

∂e 3n

, e 1 , e 2 , . . . e 3n



= U (e 1 , e 2 , . . . e 3n ) (I.18)

as H = F − U and F can be expressed in terms of the η i and the partial derivatives of S. (We discuss the second equation below.) The first of this pair, (I.17), is the Hamilton-Jacobi equation, although today it is written using H, not F and U , that is,

∂S

∂t + H  ∂S

∂η i

, η i



= 0. (I.19)

A Comment on Constraints

Hamilton writes in the Second Essay

“... it is useful to express previously the 3n rectangular coordinates x y z as functions of 3n other and more general marks of position η 1 η 2 . . . η 3n ..”

and thus explicitly says that there are as many η i as there are Cartesian co- ordinates. However it is not needed in his arguments that the number of η i

is 3n, and the reasoning still holds when there are fewer η i , as is the case for constrained systems. In the Introduction to the General Method [9], Hamilton writes that his method is applicable to

“... any system ... of attracting or repelling points, (even if we suppose those points restricted by any conditions of connexion con- sistent with the law of living force)”

so it is not the case that he thought this assumption was necessary.

Just as Hamilton believed conservation of energy would be found to be true if a seemingly non-conservative system was analysed fully, he wrote in the same introduction

“the science of force ... [tends] more and more to resolve all connex- ions into repulsion and attraction of points”

which is the modern view. (For example, [1] considers constraints as the limit of having a very steep potential around a curve or surface in space).

Hamilton’s Second Partial Differential Equation

In both the General Method and Second Essay, Hamilton presents two PDE for the characteristic function V and the principal function S, respectively. One of these equations is concered with the endpoints of the motion and the other – (I.18), which appears as (C) in the Second Essay – with the inital points.

Indeed, Jacobi criticized Hamilton for giving two equations for the same function S, partly for the lack of proof that a simultaneous solution exists, partly for the superflousness of a second equation when a solution to the first equation solves the mechanical problem [11]. Jacobi also argues that for a non-conservative system, the second equation cannot be expected to hold, but the first remains valid also for time-dependent potentials U .

In [14] the authors note that while in the “final” equation (I.17) the endpoints

x _i are functions of time, in the initial equation (I.18) time is fixed t = 0 and

the a i are not functions of time. According to them, this is why Hamilton does

not actually use the second equation in either paper. It is possible to reply to

the criticism of Jacobi by allowing to the initial time to be a variable τ as in [14] or t 0 as in [3] to make the initial equation mathematically more meaningful.

However, it is not clear what, if any, physical information is gained and it seems that a clearer physical pictures is given by one equation, the solution of which determines the motion, given the initial state.

Only the final equation (as a canonical transformation), is discussed in Arnold [1], Goldstein [7] and Landau [13]; it seems that history has come to agree with Jacobi’s assessment of the initial equation as redundant.

The Optical Analogy

The principle of least action is not the only variational principle in physics, and not the first. As indicated by the title of his report On the Application to Dynamics of a General Mathematical Method Previously Applied to Optics [8] Hamilton’s ideas originated in his works on optics. There he used Fermat’s principle: light travels between points along the path that takes the shortest time.

Using this variational principle Hamilton derives that an optical system can be described by a characteristic function V that should satisfy

X  ∂V

∂x i

 2

= n ²

where n is the index of refraction. From the derivatives of V the paths of rays can then be found. This is essentially the same as equation (I.7) since n is allowed to vary in space.

In this way we have a description of classical particle mechanics analogous to the description of classical waves. At least in notation, we have the formulation of classical mechanics closest to quantum mechanics. A notational sleight of hand and throwing in some constants gives

i~ ∂S

∂t + ˆ H

 i~ ∂

∂q , ˆ q



S = 0 (I.20)

the Schr¨odinger equation for S. Of course, the interpretation of the function S and how to extract physical information from it is different, and ˆ H a linear operator on a Hilbert space, as opposed to the (non-linear) function H on 2n- dimensional real space.

I.2 The Principal Function as a Canonical Trans- formation

A modern account of the Hamilton-Jacobi equation, as presented in for example [1, 7, 13] starts from Hamilton’s canonical equations

˙q = ∂H

∂p ˙p = − ∂H

∂q

where we call the q:s coordinates and the p:s momenta. One then considers a

transformation to a new set of coordinates Q and momenta P such that the

canonical equations take the same form Q ˙ = ∂K

∂P P ˙ = − ∂K

∂Q

where K is a new Hamiltonian, expressed in the new coordinates and momenta.

These transformations are called canonical. For instance, one could use polar or spherical coordinates instead of Cartesian, angular momenta instead of linear, or some other coordinate system.

We can formulate Hamilton’s modified principle, that the variation of the action vanish, in these two systems.

δ Z t

t

P _i · ˙Q i − K(Q, P, t) dt = 0 (I.21)

δ Z t

t

p · ˙q ⁱ − H(q, p, t) dt = 0 (I.22) Here the variations are taken to vanish at the end-points. These equations will be satisfied if we have a relation such as

p · ˙q − H = P · ˙Q − K + dF dt

where F is any function of phase space and time. ² This function F is called the generating function. We can for example take F to be a function of Q and q and time. Then it must hold that

p · ˙q − H = P · ˙Q − K + ∂F

∂t + ∂F

∂q ˙q + ∂F

∂Q · ˙ Q and since all the Q, q are independent

P = − ∂F

∂Q p = ∂F

∂q and

K = H + ∂F

∂t . (I.23)

It is also possible to let F be a function of some other combination of the new and old coordinates and momenta and similar relations in terms of the partial derivatives of F will be obtained, as listed in [7] and Table I.1. Note that not every canonical transformation can be generated with an F of one of these forms.

If we could find a canonical transformation such that transformed Hamilto- nian vanishes the equations of motion are trivial, for then both the coordinates and momenta are constant. A solution in terms of the old variables is then obtained by inverting the transformation. The condition, using (I.23) is

H(q, p, t) + ∂F

∂t = K = 0

2

There could be a scaling factor, but a rescaling of the coordinates and momenta is a

canonical transformation, so we can assume the scaling factor to be 1 [7].

Table I.1: Canonical transformations generated by generating functions.

Generating function Coordinate relations F = F 1 (q, Q, t) p = ^∂F _∂q

P = − ^∂F _∂Q F = F 2 (q, P, t) p = ^∂F _∂q

Q = ^∂F _∂P

F = F 3 (p, Q, t) q = − ^∂F _∂p

P = − ^∂F _∂Q

F = F 4 (p, P, t) q = − ^∂F ∂p

Q = ^∂F _∂P

and if we take the function F = F (q, P, t), then p = ^∂F _∂q and we obtain the Hamilton-Jacobi equation

H

 q , ∂F

∂q , t

 + ∂F

∂t = 0 (I.24)

which with a relabeling is the same as Hamilton’s equation for the prinicpal function S, (I.17).

A solution to (I.24) can have n independent constants of integration α i , i = 1, . . . , n. Since the equation involves only derivatives of S, we can add another constant, S ⁰ = S + α n+1 but this is superfluous. Such a family is called a complete integral or solution. We could write the solution as S = S(q, α i , t) and pick α i = P i , where P i are the new, constant, momenta. Since S is a canonical transformation we have the relations

p = ∂S(q, α i , t)

∂q Q i = ∂S(q, α i , t)

∂α i

.

Since we know q, p at t = 0, the first relation gives a system of equations for the α i and inverting the second relation yields the q as functions of time, solving the mechanical problem.

Suppose the Hamiltonian is conserved. Then if S solves the Hamilton-Jacobi equation, ^∂S _∂t = H = E, a constant. Then S must have the form

S(q, t, α i ) = W (q, α i ) − Et (I.25) where i = 1, . . . , n − 1 and W solves

H

 q, ∂W

∂q



= E.

This function W is called the characteristic function, and we see that it corre- sponds to Hamilton’s characteristic function V , from (I.9).

Finally we note that this view is in contrast to Hamilton’s original work, where he, quoting [14]

“assumed that the principal function S has a specific form: it is a function of time t, the initial coordinate values e i , and the final coordinate values η i . The possibility of other forms for S was not envisaged by him.”

In fact, calling momenta coordinates and vice versa is a canonical transforma-

tion, even if it is only applied to some coordinates and momenta. Goldstein

writes [7]

“There is no longer present in the theory any lingering remnant of the concept of q i as a spatial coordinate and p i as a mass times a velocity.”

This is an abstraction far from Hamilton’s original work.

I.3 Hamilton’s Equations as Characteristics

It is possible to express the relation between the system of ODE, Hamilton’s canonical equations, and the Hamilton-Jacobi equation PDE also “the other way around”. Suppose that we have a first-order PDE for u = u(x),

F  ∂u

∂x , u, x



= 0 (I.26)

in some region of R ⁿ . Then perhaps it is possible to find a trajectory x(s), parametrized by s, along which we can find u. Let

z(s) = u(x(s)) p = ∂u

∂x (x(s)).

In our context of mechanics, z will correspond to the action and p the momenta, but this method is general.

To find the path x(s) we differentiate p:

dp i

ds =

n

X

j=1

∂ ² u

∂x i ∂x j

(x(s)) dx j

ds (s)

and to eliminate second derivatives, we use that u satisfies a PDE, so that

0 =

n

X

j=1

∂F

∂p j

∂ ² u

∂x j ∂x i

+ ∂F

∂z

∂u

∂x i

+ ∂F

∂x i

.

Evidently, if

dx j

ds (s) = ∂F

∂p j (I.27)

holds, then

dp i

ds = − ∂F

∂x i − ∂F

∂z p i (s) (I.28)

An equation for z is given by differentiating the definition of z and using the previous equations. Thus

dz ds =

n

X

j=1

∂u

∂x j

dx j

ds =

n

X

j=1

p j ∂F

∂p j . (I.29)

Collecting these equations in vector form we obtain the characteristic ODE:s for the equation (I.26),



 

 

dp

ds = − ^∂F ∂x − ^∂F ∂z p

dx ds = ^∂F _∂p

dz

ds = ^∂F _∂p · p

(I.30)

Now consider the Hamilton-Jacobi equation H

 q, ∂S

∂q

 + ∂S

∂t = 0.

If x = (q, t), the n + 1:th component of the equation for ^dx _ds is ^dx _ds

= 1, since the n + 1:th component of p is ^∂S _∂t . We can take s = t and we have the characteristic equations



 

 

dq

dt = ^∂H _∂p

dp

dt = − ^∂H _∂q

dz

dt = ^dq _dt · p − H

(I.31)

which are Hamilton’s canonical equations together with z as the action of the

system.

Chapter II

Classical Perturbation Theory Applied to

Planetary Orbits

Much of Hamilton’s aim with the development of the technique of the principal function was to find methods of approximation. Mechanical problems with ex- act, simple solutions are few. However, often a more complicated problem differs only a little from one that has an exact solution. We call this a perturbation, and we want to express the modified motion in terms of the perturbation and the exact solution.

II.1 Canonical Perturbation Theory

Let H(q, p, t) = H 0 (q, p, t) + ∆H(q, p, t) where we imagine that ∆H  H 0 . If the Hamilton-Jacobi equation can be solved for H 0 , we have a function S that generates a canonical transformation to constant momenta and coordinates (α, β). This is still a canonical transformation for the Hamiltonian H, so

K(α, β, t) = H 0 + ∆H + ∂S

∂t = ∆H(α, β, t)

since S solves the Hamilton-Jacobi equation (I.24) for H 0 (equivalently, the transformation generated by S makes the transformed Hamiltonian vanish iden- tically).

Then the transformed canonical equations are

˙

α i = − ∂∆H

∂β i

β ˙ i = ∂∆H

∂α i

and these are exact. However, if ∆H is small a first order approximation is to take

˙α 1i = − ∂∆H(α, β, t)

∂β i

    ₀

β ˙ 1i = ∂∆H(α, β, t)

∂α i

    ₀

where we evaluate the right-hand sides for the constant (α 0i , β 0i ). If these

equations can be solved, the approximation can be refined by using the first-

order approximation as the starting point and repeating the process.

For example, consider a harmonic oscillator potential as a perturbation on free motion,

H = p ² 2m + k x ²

2 .

Then since the Hamiltonian is conserved, from (I.25) S = W − Et and

 ∂W

∂x

 2

= 2mE which gives

S = xα + α ² t 2m

where α is the momentum. Then the transformed constant coordinate Q is Q = ∂S

∂α = x − αt m which is clearly linear motion and we can let Q = β.

The transformed Hamiltonian is

∆H = k 2

 β + αt

m

 2

with canonical equations

˙α = − ∂H

∂β = −k(β + αt

m ) β = ˙ ∂H

∂β = kt

m (β + αt

m ). (II.1) This system can be exactly solved (and the solution is sinusoidal motion). If we instead take the perturbation approach and let the initial position be x = β 0 = 0 and initial momentum be p 0 = α 0 we have

˙α 1 = − α 0 kt

m β ˙ 1 = t ² k m ² α 0 . We can see that, using the inital conditions,

α 1 = α 0 − kα 0 t ²

2m β 1 = kα 0 t ³ 3m ² . We can put these solutions into (II.1) to obtain

˙α 2 = −k  α 0 t

m − kα 0 t ³ 6m ²



with solution

α 2 = α 0 − kα 0 t ²

2m + k ² α 0 t ⁴ 24m ²

which, with ω ² = k/m, is α 0 cos(ωt) to fourth order. One can check that β 2 is

proportional to sin(ωt) to fifth order, as expected.

II.2 Separability and Action-Angle Variables

Suppose that S is a solution to the Hamilton-Jacobi equation and S can be written in the form

S(q, α i , t) = S 1 (q 1 , α i ) + S ⁰ (q 2 , . . . , q n , α i , t) i = 1, . . . , n.

That is we can separate S into one part that depends only on q 1 , and one that depends only on the other variables. If also a solution of this form separates the Hamilton-Jacobi equation into one equation for S ⁰ and one for S 1 , we say that q 1 is a separable coordinate. Now if the Hamiltonian also does not depend explicitly on time (is conserved), we know that (I.25) S(q, α i , t) = W (q, α i )−Et since − ^∂S _∂t = H = E, which is a constant. But then, by hypothesis,

E = X

i

H i (q i , α k ) = H j (q j , α k ) + X

i6=j

H i (q i , α k )

and subtracting H j from both sides, E − H j (q j , α k ), which depends only on q j , is equal to something that does not depend on q j and must be a constant. Then each H j solves an equation

H i

 q i , ∂W

∂q i

, α j



= α i

which is a first order ordinary differential equation, so it can be reduced to quadratures. If all coordinates are separable, we say that the system is com- pletely separable.

Suppose we have a completely separable, conservative system with coordi- nates and momenta (q i , p i ). If for every coordinate-momentum pair, projecting the motion of the system in time on the (q i , p i ) plane in the phase space gives a closed curve or a periodic function in q i we can construct action-angle vari- ables. We call them so because complete separability implies that there are n independent constants of motion (one of which is the Hamiltonian). If this is the case, the region of phase space defined by fixing each constant of motion to some value (giving initial conditions) is diffemorphic to an n-dimensional torus [1]. A natural choice of coordinates on a torus is of course a set of angles, and the canonically conjugate variables will have the dimension of action.

The overall motion need not be periodic. For example, a 2-dimensional harmonic oscillator with frequency 1 Hz in the x-direction and √

2 Hz in the y-direction will trace out ellipses in both the (x, p x ) and (y, p y ) planes, but as

√ 2 is not a rational number no integer number of x periods is also an integer number of y periods.

We define the action variable by J i =

I

p i dq i (II.2)

where the integral is taken over a period of the (q i , p i ) plane orbit. It is clear

that J i has units of action. If q i is cyclic and p i constant, we take q i to go from

0 to 2π, as is natural for an angle. That is, for cyclic p i , J i = 2πp i . Note that

given inital conditions, the J i are constants, and we can take the n integration

constants α i to be functions of the J i .

W generates a canonical transformation where the Hamiltonian is cyclic in all coordinates. Remembering that differentiating with respect to the integration constants gives canonical coordinates, we have

w i = ∂W

∂J i

˙ w i = ∂H

∂J i

where H is a function only of the J j , which are constants. Hence right-hand side of Hamilton’s canonical equation for w i is always taken in the same point, that is, is constant, and the conjugate angle variables w i are linear functions of time

w i = v i t + β i

where v i is the frequency of the motion in the i-coordinate.

Note that with this definition, used in [7], after one period in the i-coordinate, w i increases by unity. We can write (for one degree of freedom)

∆w = I ∂w

∂q dq =

I ∂ ² W

∂q∂J dq = ∂

∂J I ∂W

∂q dq = 1

where W is the principal function, so that ^∂W _∂q = p, making the last integral the definiton of J. It is possible to instead define the action variable (II.2) with a factor (2π) ⁻¹ , so that the angle variables increase by 2π for each period. This definition is used by [1].

II.3 Action-Angle Variables for the Kepler Prob- lem

Consider a planet orbiting around its star due to their gravitational interaction.

The problem of determining the orbit is the Kepler problem. It is well known that bound solutions (orbits such that the planet stays in a bounded region of space) are ellipses in a plane. We will construct a solution in terms of action- angle variables. We will use units where Gm Earth = 1 and work in a center of mass system, but approximate the reduced mass with m Earth .

In spherical coordinates (r, θ, φ) the Hamiltonian is

H = T + U = 1

2m p ² _r + p ² _θ

r ² + p ² _φ r ² sin ² θ

!

− M r . Clearly φ is cyclic. A separated solution then, has the form

W = W r (r) + W θ (θ) + α φ φ.

Using this ansatz the Hamilton-Jacobi equation (I.19) is

 ∂W r

∂r

 2

+ 1 r ²

"

 ∂W θ

∂θ

 2

+ α ² _φ sin ² θ

#

− 2mM

r = 2mE

where E is the energy, which is conserved. Multiplying by r ² , we can move terms around to obtain

r ²  ∂W r

∂r

 2

− 2mMr − 2mr ² E = −  ∂W θ

∂θ

 2

− α ² _φ

sin ² θ

and as the left side depends only on r and the right side only on θ, both sides must be constants. The equation is then solved by quadratures:

 ∂W θ

∂θ

 2

= α ² _θ − α ² _φ

sin ² θ (II.3)

W θ = Z θ

θ

s

α ² _θ − α ² _φ

sin ² θ dθ (II.4)

 ∂W r

∂r

 2

= 2mE + 2mM r − α ² _θ

r ² (II.5)

W r = Z r

r

r

2mE + 2mM r − α ² _θ

r ² (II.6)

and the separation constants are identified as the z-component of angular mo- mentum (α φ ), its total magnitude (α θ ) and the total energy. We can temporarily choose the coordinate system so that α φ = α θ = l. Then W θ = 0 and one of our constant coordinates is

β = ∂S

∂l = φ + ∂W r

∂l = φ − Z r

r

l r ²

dr q

2mE + ^2mM _r − _r ^l

and this integral can be evaluated (by changing variables to u = 1/r and com- pleting the square) to find

β = φ − arccos

l

mMr − 1 q

1 + _mM ^2El

.

If we define

l ²

mM = a(1 − e ² ) E = − M 2a

(here E < 0 as we consider bound orbits) and solve for r, we find r = a(1 − e ² )

1 + cos(φ − β) (II.7)

which is the equation for an ellipse with eccentricity e and semi-major axis a.

If the angular momentum is not along the z-axis, (II.7) will still hold since the shape of the orbit cannot depend on our arbitrary choice of z-axis. The orbit will however be parametrised not by the angle φ in the xy-plane, but by an angle ϕ in the orbital plane.

Action-Angle Variables

φ is cyclic, so the action variable J φ is just 2π times the constant value of the momentum conjugate to φ, J φ = 2πα φ = 2πp φ . To find the other action variables we use p θ = ^∂W _∂θ , so

J θ = I ∂W

∂θ dθ = I

s

α ² _θ − α ² _φ

sin ² θ dθ

where we use the separated Hamilton-Jacobi equation to solve for ^∂W _∂θ . This integral can be solved with trigonometric substitution and the result is [7]

J θ = 2πα θ (1 − cos i) = 2π(α θ − α φ ) where cos i = ^α _α

, the angle of inclination of the orbit.

The integral for J r is J r =

I r

2mE + 2mM r − α ² _θ

r ² dr which is more difficult to solve [7], but gives

J r = −(J θ + J φ ) + πM r 2m

−E (note that E < 0 for bound orbits).

For bound orbits all frequencies are equal since the orbit is a closed curve (an ellipse). This can be confirmed by writing the Hamiltonian as

H = E = − 2π ² mM ² (J r + J θ + J φ ) ²

and seeing that it depends only on the sum of all three action variables, so that the canonical equations all have the same form:

˙

w r = ∂H

∂J r

= ∂H

∂J θ

= ˙ w θ = ∂H

∂J φ

= ˙ w φ .

To summarize, we have found a solution S to the Hamilton-Jacobi equation for the Kepler problem: S = W − Et, where

W = φα φ + Z θ

θ

s

α ² _θ − α ² _φ sin ² θ dθ +

Z r r

r

2mE + 2mM r − α ² _θ

r ² and the action variables in spherical coordinates

J φ = 2πα φ J θ = 2π(α θ − α ^φ ) J r = −2πα ^θ + πM r 2m

−E .

Orbital Elements in Terms of Action-Angle Variables

It is easier to relate certain linear combinations of the angle variables w φ , w θ , w r

to the geometry of the orbit. Using the linear transformation w 1 = w φ − w ^θ

w 2 = w θ − w r

w 3 = w r

and the conjugate action variables J 1 = J φ

J 2 = J φ + J θ

J 3 = J φ + J θ + J r

Table II.1: Angles used in the text. See also Figure II.1 and Figure II.2.

Symbol Meaning

α = φ The angle relative to the x-axis of the projection of Earth’s position vector on the xy-plane ψ The eccentric anomaly of the Earth ϕ The true anomaly of the Earth

i The inclination

Ω Longitude of the ascending node. Proportional to the canonical angle variable w 1 .

ω Argument of the perihelion. Proportional to the canonical angle variable w 2 .

β The angle between Jupiter’s position vector and the x-axis (the true anomaly)

η The eccentric anomaly of Jupiter

w 3 The mean anomaly of Earth. Canonical angle variable.

W 3 The mean anomaly of Jupiter. Canonical angle variable.

we can see that the frequencies of w 1 and w 2 are zero. The constant values of these angle variables will describe the orientation of the orbit. We will relate w 3 to the position along the orbit in the section below.

Looking at Figure II.1, ω = 2πw 2 is the argument of the perihelion [7].

Ω = 2πw 1 is the angle to the line of nodes. J 1 = J φ is the z-component of angular momentum. J 2 = 2πl where l is the total angular momentum. The semi-major axis a is given by

a = J ₃ ² 4π ² mk and is independent of J 1 and J 2 .

Using the angles as defined in Table II.1 (see also Figure II.3), the following relations hold [7].

2πw 1 = Ω (II.8)

2πw 1 + u = α (II.9)

2πw 2 = ω (II.10)

sin u = cot i cot θ (II.11)

cos θ = sin i sin(ϕ + ω) (II.12)

Note that the angles are defined in terms of their sines or cosines and we will

typically only need trigonometric functions of them. It is then not necessary to

compute the actual value of u, for cos u and so on can be found with trigono-

metric identities. This will, however, often yield unwieldy expressions and we

will not always express everything in the “fundamental” angles (Ω, i, ω, ϕ) that

give the orbit and the position along it.

x

y z

Ω ω i

ˆ n

Figure II.1: Orbital elements. The vector ˆ n is normal to the orbital plane.

x

y z

N P

R

α

Ω u

θ

Figure II.2: Angles used for the Kepler problem. ON is the line of nodes; R is

the position of the orbiting body.

True, Eccentric and Mean Anomalies

The Kepler problem has three degrees of freedom, so we need six numbers to describe the state of the system. w 1 , w 2 and i = cos ^{−1 J} _J

2

determine the orientation of the ellipse. The semi-major axis and eccentricity of the ellipse can be determined in terms of J 2 and J 3 . We need finally a number to tell the position along the ellipse and this must involve w 3 since we haven’t used it for any of the other five numbers.

A natural choice to tell the position along the ellipse is the angle relative to the axis between the periapsis and the focus of the ellipse, where the sun is, as in Figure II.3. This is the true anomaly ϕ. But w 3 is a canonical angle variable so it is traversed with constant frequency and from the conservation of angular momentum l = mr ² ϕ this is only the case for ϕ for a circular orbit. So ˙ in general ϕ that most easily describes the motion in real space is not a simple function of w 3 , that most simply describes the motion in phase space. This will be the source of complications.

F

ψ ϕ

O

Figure II.3: Eccentric anomaly ψ and true anomaly ϕ.

Let a and b be the axes of the ellipse. Draw a circle C with radius a concentric with the orbital ellipse. Let ψ be the angle from the center to the point on C that lies directly above the position of the planet, as in Figure II.3. We call this the eccentric anomaly ψ. The eccentricity is e =

q

1 − ^b a

and if the origin

is at the center of the ellipse, the foci are at (±ae, 0). If we for the moment

let x, y be coordinates in the orbital plane, putting the origin at the center

of the ellipse, by definition x = a cos ψ, and from the equation for an ellipse,

x ² /a ² + y ² /b ² = 1 we conclude y = b sin ψ. Then by the Pythagorean theorem

r ² = b ² sin ² ψ + (a cos(ψ) − ae) ² . That is,

r ² = b ² sin ² ψ + a ² e ² + a ² cos ² ψ − 2a ² e cos ψ

= b ² sin ² ψ + a ² − b ² + a ² cos ² ψ − 2a ² e cos ψ

= a ² − b ² cos ² ψ + a ² cos ² ψ − 2a ² e cos ψ

= a ² + a ² e ² cos ² ψ − 2a ² e cos ψ = a ² (1 − e cos ψ) ² (II.13) If we take E to be a momentum, the conjugate coordinate is, from the equations of transformation and (II.6),

Q E = ∂W

∂E = ∂W r

∂E = r m 2

Z r r

dr q

E + ^M _r − _2mr ^l

but the canonical equation here is ˙ Q E = ^∂H _∂E = 1, so if we take the starting time to be when r = r 0 , and use the relations _mM ^l

= a(1 − e ² ) and E = − ^M 2a for elliptical orbits [7] we can write

t =

r m

2M Z r

r

r dr q

r − ^r _2a

− ^a(1−e ₂

⁾

and change variables to ψ, using (II.13), to get after some algebra

t = r ma ³

M Z ψ

0 (1 − e cos s) ds = r ma ³

M (ψ − e sin ψ)

so the quantity ψ − e sin ψ is a linear function of time and goes from 0 to 2π in one period. We conclude that

2πw 3 = v 3 t = r ma ³

M t = ψ − e sin ψ (II.14) which is called Kepler’s equation and 2πw 3 is called the mean anomaly.

The orbit equation (II.7) can then be written in two forms:

r =

( _a(1−e

₎

1+e cos ϕ = _mM ^l

_{1+e cos ϕ} ¹

a(1 − e cos ψ) = _mM ^l

^{1−e cos ψ} _1−e

. (II.15) where ϕ is the true anomaly, ψ is the eccentric anomaly. Comparing the two forms we can find that the eccentric anomaly ψ is related to the true anomaly ϕ by

cos ϕ = cos ψ − e

1 − e cos ψ . (II.16)

II.4 Perturbation

Consider a second planet which we take to be much more massive than the first,

but much less massive than the star. We shall call the planets Earth and Jupiter

to reflect the hierarchy of masses. This hierarchy of masses means that Jupiter

affects Earth’s orbit much less than the sun does, and Jupiter is approximately

uninfluenced by Earth. We thus consider Jupiter as a perturbation, orbiting in an ellipse. Its orbit can thus be described with the same set of action-angle variables. We will only need the angle variable W 3 , the mean anomaly.

The Hamiltonian for the perturbed system is H = T − M

r − M ⁰

d (II.17)

where T is the kinetic energy and d is the distance between Earth and Jupiter, which has mass M ⁰ . The last term is small; it is our perturbation ∆H. Let r(ψ), R(β) be the orbits of the planets. Then

d ² = r ² + R ² − 2rR cos(α − β) sin θ (II.18) where α is the angle in the xy-plane of Earth and β the same for Jupiter, which is its true anomaly, and θ is the polar angle in spherical coordinates (0 ≤ θ ≤ π/2).

II.5 Averaging

We seek the average of the perturbation of the perihelion (or another orbital parameter). We could take the average over one period of the motion

˙ω = 1 τ

Z τ 0

∂H

∂l dt

(where we used 2πw 2 = ω, and 2πJ 2 = l) but ^∂H _∂l is a function of both positions and the unperturbed positions are both periodic, with different periods, ¹ so it is not clear which value should be used for τ . We would like to average over all possible configurations, so we can take the limit when τ → ∞. Then the time average can be replaced with a space average [1] in the sense that we average over the space of canonical angle variables (a torus). Thus (where angle variables go from 0 to 1)

˙ω = lim

τ →∞ τ ⁻¹ Z τ

0

∂H

∂l dt = Z 1

0

Z 1 0

∂H

∂l dw 3 dW 3 (II.19) But from Kepler’s equation (II.14) we have

∂(w 3 , W 3 )

∂(ψ, η)    

= (2π) ⁻²    

1 − e 1 cos ψ 0 0 1 − e 2 cos η

=

= (2π) ⁻² (1 − e ¹ cos ψ)(1 − e ² cos η)

where ψ, η are the eccentric anomalies of Earth and Jupiter, respectively. Thus we can find the average perturbation by integrating over the eccentric anomalies, in terms of which we can express trigonometric functions of other angles.

Note that if we had used angle variables that go from 0 to 2π, the factor (2π) ⁻² would have shown up in (II.19) instead of in the functional determinant.

1

This is not a problem if the periods are rationally commensurate τ

/τ

∈ Q for then we can choose τ as the least common multiple and the whole system is periodic with this period.

In fact what follows assumes that the frequencies are not rationally commensurate.

II.6 Perturbations of Orbital Elements

Perturbation of the Perihelion

Consider the perihelion ω = 2πw 2 . The action variable conjugate to w 2 is J 2 = 2πl, so the perturbation of ω is

˙ω = ∂∆H

∂l = ∂

∂l

 −M ⁰ d



which we can compute using the orbit equations for the planets. Explicitly, we have

∂

∂l 1 d = ∂

∂l l ⁴ m ² M ²

 1 − e cos ψ 1 − e ²

 2

+ R ² − 2R l ² mM

1 − e cos ψ

1 − e ² cos(α − β) sin θ

! −

where we have not written out R since it is independent of l. The derivative is

∂

∂l 1

d = − 4 l ³ m ² M ²

 1 − e cos ψ 1 − e ²

 2

− 4 l

mM R 1 − e cos ψ

1 − e ² cos(α − β) sin θ

! 1 2d ³

= − 2

l (r ² − Rr cos(α − β) sin θ)d ⁻³ .

To take the average, we use relations (II.8)–(II.12) to express α, β in terms of ψ, η, and average over ψ and η. That is,

˙ω = 1 (2π) ²

2M ⁰ l

Z Z r ² − Rr cos(α − β) sin θ

d ³ (1 − e 1 cos ψ)(1 − e 2 cos η) dηdψ (II.20) where the integration is over [0, 2π] × [0, 2π].

The integration can be done numerically and using values for the Sun, Earth and Jupiter [15] the average precession rate is around 20 seconds of arc per year. With a computer we can easily vary the paramaters and in Figure II.4 we show values for varying the eccentricity of the orbit of Jupiter between 0.05 and 0.5. The average precession rate decreases with increasing eccentricity.

Physically we can reason that as the eccentricity increases, the orbit is flattened and from Kepler’s second law, or the law of areas, or the conservation of angular momentum, the planet moves more slowly when it is further away from the star.

This means that more time is spent where the perturbing force is smaller.

It is also possible to solve the (Newtonian) equations of motion numerically

with a computer. Figure II.5 shows this for coplanar orbits. To present the effect

(which realistically is very small) we have set the perturbing planet’s mass to

half that of the star rather than a difference of three orders of magnitude. The

precession of the periapsis is then seen even after a few orbits. Note that since

the perturbation is exaggerated, the shape of the orbit is distorted and does not

obviously resemble a rotating ellipse. We should also note that with the values

used here, the perturbing planet, orbiting in a circle with radius 50 units, orbits

much more slowly than the smaller planet. In Figure II.5 the perturbing planet

has only moved slightly less than a quarter around the star. Thus the perturbing

force lies in the first quadrant for the plotted time.

0 0.1 0.2 0.3 0.4 0.5 0

2 4 6 8 10 12 14 16 18 20 22

Eccentricity

Average precession rate [arc seconds per year]

Precession of the perihelion

Figure II.4: Average precession of the perihelion using values for the Sun, Earth and Jupiter, varying the eccentricity of Jupiter’s orbit.

−5 0 5 10

−6

−4

−2 0 2 4 6 8 10

x

y

Orbit 1 Orbit 2 Orbit 3 Orbit 4 Orbit 5

Figure II.5: Trajectory for a planet perturbed by another planet with half the

mass of the star, orbiting in a circle with radius 50 (units arbitrary). The star

is fixed at the origin (yellow dot). Initially, both planets are on the x-axis and

orbit counter-clockwise. To aid visualization different colors are used along the

orbit.

Perturbation of the Ascending Node

We have that 2πw 1 = Ω. The canonical variable conjugate to w 1 is J 1 , so we should have

˙Ω = 1 2π

∂∆H

∂J 1

.

But J 1 = J 2 cos i, so the precession in Ω is simply proportional to the precession in ω:

˙Ω = 1 cos i ˙ω which is geometrically reasonable.

Perturbation of the Inclination

The inclination i is given by cos i = ^J _J

, so

−˙i sin i = J ˙ 1

J 2 − J 1 J ˙ 2

J ₂ ² = J ˙ 1 − cos i ˙ J 2

J 2

and we can assume that sin i 6= 0, since if this were not the case the motion would be planar and ˙i = 0. Of course,

J ˙ 1 = − ∂H

∂w 1 = −2π ∂H

∂Ω J ˙ 2 = − ∂H

∂w 2 = −2π ∂H

∂ω and we can realize that

∂

∂ω 1 d = 2

d ³ Rr ∂(cos(α − β) sin θ)

∂ω since r and R are functions only of ψ and η. Likewise

∂

∂Ω 1 d = 2

d ³ Rr ∂(cos(α − β) sin θ)

∂Ω .

Now, from relations (II.8)–(II.12), we have ²

∂α

∂ω = ∂u

∂ω cos u ∂u

∂ω = cot i 1 sin ² θ

∂θ

∂ω

− sin θ ∂θ

∂ω = sin i cos(ϕ + ω) and we also have ^∂α _∂Ω = 1. This gives

∂α

∂ω = cos i

cos u sin ³ θ cos(ϕ + ω).

Then

∂

∂w 1 d = 2Rr

d ³



cos(α − β) cos θ sin i cos(ϕ + ω)

sin θ − sin θ sin(α − β) cos i

cos u sin ³ θ cos(ϕ + ω)



2

We take H = H(J

, J

, J

, w

, w

, w

). Then i = i(J

, J

), ψ = ψ(w

) and u =

u(i, w

, w

) so that, in particular,

= 0.

= 2Rr

d ³ cos(ϕ + ω)



cos(α − β) cot θ sin i − sin(α − β) cos i cos u sin ² θ

 . Now similarly

∂

∂Ω 1 d = 2Rr

d ³



− sin(α − β) sin θ dα

dΩ + cos(α − β) cos θ dθ dΩ

 .

We can realize that since Ω corresponds to a rotation around the z-axis, changing Ω has no effect on the polar angle θ, that is _∂Ω ^∂θ = 0.

We can use relations (II.8)–(II.12) to simplify this somewhat. Note that cos i sin i cot θ = cot θ cot i sin ² i = sin u sin ² i.

A full expression for the perturbation of the inclination is then

− di

dt sin i = 4πM ⁰ Rr J 2 d ³



− cos(ϕ + ω)



cos(α − β) sin u sin ² i−

sin(α − β) cos ² i cos u sin ² θ



− sin(α − β) sin θ

 (II.21)

and since i, and thus sin i are typically very small we might worry that, dividing by sin i, ^di _dt would blow up. However looking into it in detail, for small i there are no zeroth or first-order terms on the right-hand side. We present the details in Appendix A.

II.7 Tidal Locking

In the Kepler problem the orbiting body is treated as a point particle. Of course celestial objects such as moons are not point particles and different parts of the body are subject to different forces. If there is sufficient asymmetry in the body there will be a net torque. For instance if two masses m and 2m are fixed at the ends of a rod initally perpendicular to the orbit the force on the heavier end will be larger and there will be a net torque.

In this section we assume a circular orbit.

Let L be the volume occupied by the body (the moon) and M the mass of the central body (the planet). Its potential energy is

U = − Z

L

GM ρ R dL

where ρ is the density and R the distance from a point in L to the main body (which we treat as a point particle or at least spherical). Now this can be expanded in terms of Legendre polynomials, where α is the angle relative to the line between the centers of mass, as in Figure II.6:

U = − GM m

r − GM

r ² Z

ρr ⁰ cos α dL + GM 2r ³

Z

ρr ⁰² (1 − 3 cos ² α) dL + . . . . The first order term vanishes if the origin is chosen as the center of mass of the orbiting body since rr ⁰ cos α = r · r ⁰ and integrating ρr ⁰ gives the center of mass. ³

3

This is the same as there being no true dipoles in gravitation.

r ⁰

α _b

r

Earth Moon

Figure II.6: An orbiting body of finite extent. (Not to scale.)

Now 1−3 cos ² α = 3 sin ² α−2 and in terms of the inertia tensor 2 R r ⁰² ρ dL = Tr I = I 1 + I 2 + I 3 where I i are the principal moments of inertia. Now r ⁰² sin α is the distance to the centers of mass axis, so this term contributes the moment of inertia around that axis. That is,

U = − GM m

r + GM

2r ³ (3I r − I 1 − I 2 − I 3 ).

Let us assume that the orbiting body is a spheroid with I 3 corresponding to the symmetry axis. Then the moment of inertia I r is I 1 + (I 3 − I ¹ )γ ² where γ = ˆ I 3 · ˆr [7]. ⁴

We need some suitable set of angles to describe the orientation of the moon.

We shall use Euler angles (φ, θ, ψ) corresponding to rotation by φ around the axis of symmetry, then rotation by θ around the I 1 axis, then ψ around the axis of symmetry. These actually correspond to the orbital elements (Ω, i, ω) as in Fig- ure II.1. Then the I 3 axis in space coordinates is (sin φ sin θ, − cos φ sin θ, cos θ) and the r-axis is (cos η, sin η, 0) where η = ωt gives the position along the orbit, so from taking the dot product γ = sin θ sin(φ − η) .

Now if φ = η + v where v is a constant the symmetry axis will precess with the orbital frequency. If also ˙ ψ = 0, the same face of the moon will be seen from the planet at all times.

By writing the angular velocities around the principal axes in terms of ˙ φ, ˙θ, ˙ ψ the kinetic energy is

T = I 1

2 ( ˙θ ² + ˙ φ ² sin ² θ) + I 3

2 ( ˙ ψ + ˙ φ cos θ) ²

and since the potential doesn’t depend on the angular velocities the conjugate momenta are

p ψ = ∂L

∂ ˙ ψ = ∂T

∂ ˙ ψ = I 3 ( ˙ ψ + ˙ φ cos θ) p θ = I 1 ˙θ

p φ = I 1 φ sin ˙ ² θ + I 3 ( ˙ ψ + ˙ φ cos θ) cos θ = I 1 φ sin ˙ ² θ + cos θp ψ

4

We can realize this by transforming to a system where the x

-axis is along α ˆ I

+ β ˆ I

+ γ ˆ I

, where ˆ I

are the principal axes, and looking at the 11-element of I in this system: I

= A

I

A

= I

A

A

= P

I

A

where A is the transformation matrix from unprimed

(principal) to primed system. But the x

axis is (α, β, γ) in the primed system, so A

=

A

= (α, β, γ).