Gravitational Potential and Tidal Forces

We can write a general expression for the gravitational potential, U, due to a body of arbitrary shape and density distribution by treating it as an accumulation of a large number of small masses, dm, each contributing potential dU, and then integrating over the masses.

           U = GM/r

This same equation can be used for a uniform sphere of mass M, where r is taken as a point outside the sphere. (i.e. the effect of the sphere is the same as if the mass were concentrated at its centre). For a point at a radius r inside the sphere, the mass M to be used in the equation is that if the material inside the radius r (parts of the sphere outside r have no effect - notice the analogy with electrostatics theory).

The same equation above can be used for the potential outside a spherically symmetric sphere which is a rough approximation for the case of the earth.

If we want to generalise, however, and look at arbitrary shapes and density distributions, we have to integrate over the effve to integrate over the effects of each small part of the body. Let us take a body of arbitrary shape in a (x,y,z) coordinate system:

The potential at a point P, distance r and direction cosines l,m,n from the centre of the coordinate system, is the sum of the effects of all the elements of the arbitrarily shaped body. One such element, P', is at position x,y,z, distance r' from P and [ro] from the centre. Thus U is given by:

We can substitute for [ro] and [psi] using the direction cosines and (x,y,z) values, and expand the bracket using the binomial expansion. Then to order [ro]²/r² we have:

So far we have kept our derivation totally general. We can take O as the centre of mass without losing generality. This is defined as:

We can also take Oz as the polar axis without losing generality. Then we have the moment of inertia about Oz = C given by:

We will now examine the simplification we can bring about by assuming a degree of symmetry in the attracting body, and at this point in the calculation we lose generality in favour of simplicity. We start by assuming Oz is an axis of symmettry. Then the momn axis of symmettry. Then the moments of inertia about Ox and Oy (A and B) will be equal:

giving:

We have also reference axes as Principal Axes:

which leads us to the simplification:

but l²+m²+n² = 1, so l²+m² = 1 - n² and so:

        U = GM/r + G/2r3[(C-A) + 3n2(A-C)]

          = GM/r + G/2r3(C-A)(1-3n2)

now,

so

where J₂ = (C-A)/MR_e²

R_e is the equatorial radius.

In general:

P₂, P₃ etc are Legendre polynomials.

GM can be obtained from the orbit of a natural satellite or an orbiting spacecraft. Any planet with a natural moon gives a method for getting the mass. For the Earth GM = 3.98603 (+/- 3) 10¹⁴ m³s^-2. In many cases for the planets the value of GM is known more accurately than M itself because of the uncertainty in the value of g n M itself because of the uncertainty in the value of g the gravitational constant, and this is the limiting factor!

The values of J_i for the earth were first determined from surface measurements. This gave J₂ and an uncertain J₄ - uncertain for the same reason that J₃ was not quoted; there was a lack of stations in the southern hemisphere (and there was insufficient data).

Since the advent of satellites there has been a 20x improvement in accuracy of the J_i's and the odd terms are now also obtained.

Thus we have

     J2 = 1082.64 (+/-0.1) 10-6     J3 =   -2.56(+/-0.07) 10-6
     J4 =   -1.52 (+/-0.3) 10-6     J5 =   -0.15(+/-0.03) 10-6
     J6 =    0.57 (+/-0.2) 10-6     J7 =   -0.44(+/-0.03) 10-6

All values of J_n (n>2) are of O(10^-6). How many J's can be determined depend on how many satellites there are and how they cover the planet.

Note that this is the relative size of the geopotential term, it is not necessarily the shape of the body - a density variation can cause these variations in the geopotential without the body being mis-shapen at all. (Gnereally it will be a mixture of shape and density variations.)

Earth Tides

We can look at the deformation of a primary body being acted upon by a secondary:

The gravitational attraction at A due to M depends on R'. The acceleration due to the moon, g_A, is (g_Ar,g_Ap) in the r and psi directions shown on the diagram. Now if P is the angle psi shown on the diagram, and P' the angle psi-dashed:

     gAr = (GM/R'2) cos P'   gAp = -(GM/R'2) sin P'

and the Tidal Acceleration is the difference between the acceleration of the point A and the acceleration of the earth's centre due to the moon. Hence:

  del-gr = gAr - gEr = (GM/R2)[(R2/R'2) cos P' - cos P]
  del-gp = gAp - gEp = (GM/R2)[sin P - (R2/R'2) sin P']

Eliminate R' and P' using RcosP = R'cosP' + r and RsinP = R'sinP' plus R'² = R² + r² - 2RrcosP

so:UP> - 2RrcosP


so:
     del-gr = (GM/R2)[R2(RcosP-r)/R'3 - cosP]
            = GM[(RcosP-r)/(R2 + r2 - 2RrcosP)3/2 - cosP/(R2]
            = GM/R3[(RcosP-r)/(1+r2/R2 - 2rcosP/R)3/2 - RcosP)
            = GM/R3[(RcosP-r)(1+3rcosP/R) - RcosP]  to O(r/R)
            = GM/R3(3rcos2P - r) to order r/R
            = GMr/R3(3cos2P-1)

     del-gp = (GM/R2)[sinP - sinP/(1 + r2/R2 - 2rcosP/R)3/2]
            = GM/R3(RsinP - [RsinP][1 + 3rcosP/R])  to order r/R
            = -3(GM/R3)rsinPcosP

So, we can say:

Where

Notice that these terms have a semi-diurnal nature. They are symmetric about the 0-180 degrees line

As the earth rotates an observer is carried +,-,+,- with respect to an Earth circular in cross-section. Note in the plane perpendicular to the diagram there is no tidal effect. Thus, though the shape is the same as that given by the J₂ gravitational term above, it is 2-dimensional rather than symmetric about the z-axis.

Chans an elastic body

Changes in g give a deformation similar to the above. (Note we have calculated the changes in potential height: this will only translate into a change in shape if the body is deformable.) If perfectly elastic the elongation axis would be towards the deforming body; otherwise there will be a phase lag.

For the solid earth the phase lag is small, though not for the sea tides. The reason for the small solid earth lag is because the period of natural oscillation is <<0.5days (it's actually around 57 minutes) so the earth can continuously and rapidly adjust. The natural period of the ocean tides is several days.

If we assume the oceans could adjust then the level of water would rise to a point of static equilibrium.

If the marine tide = del-H then the work done against gravity would be g(del-H)=-W

For the solid earth work is also done in the elastic deformation and so the height of the earth tide delH_s is less than del-H

If del-H_s = h(del-H) = h(-W/g) then h is independent of r and psi and is called the Love Number.

The liquid surface remains an equipotential. On a solid surface there is an additional potential due to the deformation of the earth and the consequent redistribution of the mass. We expectthis to be proportional to W(r,P) - say kW. k is another Love number.

A liquid surface covering the globe would remain ring the globe would remain an equipotential and be lifted (1+k)W/g relative to the centre of the earth,or (1+k-h)W/g relative to the sea bed. We can measure h and k by observing g at the surface.

      g = g0[1-(2W/rg0)(1 - 3h/2 + k)]

W is at the surface a semi-diurnal variation in g pf c. 2 in 10⁷. Changes in g of 1 in 10¹⁰ or 1 in 10¹¹ can be measured. (1 + 3h/2 + k) is c. 1.1 to 1.26%. Vertically (1 + k - h) = 0.54 - 0.82. This then gives k = 0.28 and h = 0.6. Thus the rise and fall of the earth's surface is half that of the ocean - several metres!

If we imagine a satellite in a prograde orbit about a (pro-grade) rotating planet:

A tidal bulge is set up on both bodies. If the planet rotates faster than the satellite orbits then the bulge can be carried ahead of the sub-satellite point by the primary's rotation:

(The angle L is usually small, of order 1 degree.) The bulge exerts a force on the satellite which has a tangential component speeding it up and causing it to move to a higher orbit:

The satellite in return slows the planet's spin rate:

So, the moon moves away from the planet and the planet rotation rate slows until either the satellite is lost or the spin and orbital periods become the same. Much of the energy is dissipated in the bodies as heat. Thus we see Io, the closest body to Jupiter, and to some extent Europa, the next, subject to large amounts of tidal heating. This is probably what supplies the power for the volcanic activity on Io, and may be keeping some of the ice on Europa liquid.

Whenever the satellites come near to conjunction the outer satellite will get a "tug" in the same direction on every orbit. The overall effect of these increments - which will be in one direction on one side of the orbit and in the other direction on the other side - will be toe other side - will be to stabilise the outer satellite into a resonant period with the inner satellite. (See Lewis for more details).

A satellite cannot approach its primary too closely (the Roche limit) or stray too far (the instability limit) without break-up or loss to the system respectively. The Roche Limit is the point within which a satellite would be torn apart by tidal forces if the only force hloding it together is its self-gravity. Taking a primary mass M, radius R, a distance d from its satellite mass m radius r (where M>>m) we have:

If the satellite is large enough (r greater than about 500km) its self-gravity dominates other cohesive forces. It will be torn apart if it approaches the primary closer than:

which is known as the Roche Limit, after the mathematician who first derived it. For our Moon d = 2.44(5.5/5.3)^1/3 or about 2.9R_E - around 18,500 km.

Roche's calculation was complicated by the fact he took into account the tidal distortions in the satellite just before break-up. We can get a similar (though not the "full") term if we just consider a rigid spherical satellite. Then the centripetal acceleration of the orbit is w²d where w is the angular velocity (radians/s)of the orbital motion. The gravitational accelmotion. The gravitational acceleration from the primary is GM/d², so:

     w = (GM/d3)1/2

The differential gravitational acceleration between the centre of the satellite and the outside edge is

     A = (GM/d2) - GM/(d+r)2  = approx 2GMr/d3

The differential centripetal acceleration B between the same two points is:

     B = w2(d + r) - w2d = w2r = GMr/d3

So

      A + B = 3GMr/d3

This must be balanced by the satellite's self-gravity GM/r². Disruption occurs at d = r(3M/m)^1/3.

Given the densities of the primary and satellite:

then d becomes

where as we can see we have the right form of the equation but the wrong constant because of our simplification.

The large natural satellites orbit beyond the Roche Limit of the primaries, though some of Saturn's moonlets are within it. These later must have other adhesive forces holding them together, as a low-orbiting artificial satellite of earth would. Saturn's rings are entirely within the Roche Limit, as are those of Jupiter, Neptune and Uranus.

All the gas giants have rings, though only Saturn's are opticallave rings, though only Saturn's are optically thick seen from earth. (Beware the beguiling solidity of the images of the Jovian, Neptunian and Uranian rings seen on Voyager imagery - these pictures have usually been image-enhanced to bring out the contrast. The rings are nothing like as bright or obvious, even in-situ, as they seem on these photographs.) The rings of Uranus and Neptune were discovered telescopically before space probes reached the planets, but it was only with the details Voyager images that we could appreciate their complex natures. The rings in detail are often seen to have fine detailed structure, with well-defined gaps and even in at least one case a "braided" structure. These fine features are often the result of resonance interactions with satellites either embedded within the ring system, or just outside the edges of it. These satellites are said to move in marshalling and shepherding orbits. See Lewis for images showing this behaviour in more detail.

This same interaction between rings and satellites also explains why, of example, the outer edge of the Jovian ring is sharp, while its inner edge is much more diffuse. It probably also explains the "spoke-like" structures seen in Saturn's rings. A detailed comparison of the ring systems of all the major planets (see, fig VI.23 of Lewis for full details) shows that the outer edge of Jupiter is coincident with the sub-satellite 1979J1 at about half the orbiat about half the orbital distance from the primary as the "classical" moon Amalthea. Saturn's F ring is shepherded by satellites 1980S26 to 1980S28. (It is this ring which is "braided".) Saturn's G ring is between 1980S1-S3 and Mimas, while the satellite Enceladus is embedded in the diffuse E ring. Uranus' nine rings are all inside the innermost satellite 1986U8. In the case of Neptune there are subsatellites both embedded in, and forming the edge of, its four rings.

Note that the smaller sub-satellites discovered by the Voyager spacecraft are often grouped together (like 1980S26-28) in resonant orbit groupings. It is difficult to tell if these satellites are of independent origin, or the remnants of the same bodies that broke up to form the rings. (It seems most likely that the rings formed due to break-up of bodies coming within the Roche Limit rather than their being the residue of material left in the system which failed to accrete - the rings' lifetimes are too short for the latter.)

The interesting structural differences between the rings of the different planet's are likely to exercise the minds of planetary dynamicists for some years to come. (For example, Saturn's rings are broad with small gaps between them, whereas Uranus' are like the negative image of that - thin rings separated by large gaps. Uranus' rings are also eccentric and not exactly coplanar.)

If we have a satellite mass m, at a distance d from its primary mass M₁, then a second large body of mass M₂ distance D from M₁ produces a perturbing acceleration on m of:

    A - B = [GM2/(D-d)2] - GM2/D2

  or approximately   2GM2d/D3

between orbiting body and primary (d << D). The gravitational aceleration b is given by:

     b = GM1/d2

Thus, the instability limit, when the perturbation starts to be of the same order as the primary's effect, is given by:

     d = (M1/2M2)1/3D

as long as M₁ << M₂. If M₁ >> M₂ we must use the exact relation:

    d3(2D-d) = (M1/M2)D2(D-d)2

If we take the earth-moon system with the sun as perturber then d is about 1.7 10⁶ km, four times the current orbital radius.

The equations of motion of two bodies moving under their mutual gravitational attraction are soluble algorithmically. However, as soon as a third body is added to the scenario, the equations become generally insoluble algorithmically and we are fonsoluble algorithmically and we are forced to numeric solutions. We shall look at the case of a third body of negligible mass moving in the gravitational field due to two larger bodies: this would be the case for, say, a space probe or a small satellite moving under the influence of the Earth and Sun.

We take "universal" axes X,Y,Z, in which one large body, mass m₀ is the centre of coordinate system x,y,z, centred on m₀ and parallel to X,Y,Z (see diagram). At a distance r₁ from m₀ is another large body, mass m₁. A small body of mass m, negligible compared to m₀ and m₁, is at distance r from m₀ and [ro] from m₁.

Now the equation of motion of m can be broken up to three components. In the X direction we have:

       m d2X   =   G m m0 (X0 - X)  +  Gm m1 (X1 - X)
         __         _____  ______      ____  ______
         dt2          r2      r         [ro]2   [ro]

while for m0 the equation of motion is:

       m0 d2X0   =   G m m0 (X - X0)  +  Gm0 m1 (X1 - X0)
          __         _____  ______      ____  ______
          dt2          r2      r       2      r          r12    r1

We can now divide the first by m, and the second by m₀, and subtract the second from the first, (noting that x = X - X₀, y = Y -Y₀ etc) to get:

 
         d2x  =  -G(m0 + m)x  +  Gm1(x1 - x)  -  Gm1 x1
         __       _______        _______         ____
         dt2          r3            [ro]3          r13

so

         d2x  +  G(m0 + m)x  =  Gm1(x1 - x)  -  Gm1 x1
         __       _______        _______       ____
         dt2         r3           [ro]3         r13

                              =  0  if m1=0

This last case where m₁=0 takes us back to the simple two-body problem. Where this is not equal to zero however, we have:

where

R is called the Disturbing Function. We thus now have an equation which looks like the simple two-body problem but with a disturbing factor built in which we can use in numerical solutions. In the case of a body moving through the solar system in the Sun's gavitational field, bute Sun's gavitational field, but susceptible to disturbances from other planets we can take them all into account by adding a disturbing function for each, and then we replace m₁ by m_i and sum over i. We start with initial position and velocity and then iterate forward numerically. The orbits of the planets are calculated like this, as are the comets (Cowell's method). If over a long stretch of orbit the integration is over something very near to an ellipse (i.e un-acted upon over a long stretch) then there is another method called Encke's method which is used.

Although the three-body problem is not generally soluble algorithmically, we can tell much about the behaviour in restricted situations. The usual "restricted" three-body problem assumes that the third body (the motion of which we are interested in) is small in comparison to the other bodies as in the example above, but here the two main bodies are comparable in size and we look at the motion of the third in the frame rest of the two-body system (which must be a rotating frame where the bodies rotate about their common centre of gravity).

So, taking axes centred on the centre of gravity of the system, containing bodies mass m₁ and m₂, a rotating frame in which m₁ and m₂ are fixed, we have a small body m very small compared to m₁ and m₂ and m₂, at position distance r₁ from m₁ and r₂ from m₂:

The energy equation, given that m is at (eta, epsilon) in this frame, is given by:

where the first term is the kinetic energy, the second the energy due to the centrifugal force and the next two the potential energy terms. Note the use of the rotating frame means we have to add the second term to the usual KE/PE energy equation.

There is no term due to the "Coriolis" Force as it acts at right angles to the motion and hence does no work (and does not change the overall energy). Putting the "constant" equal to:

For a given value of the parameter we can define limits to the areas m can be in given its initial energy. The defining boundary will be v²=0, where v is m's velocity. This is because v²=0 is the maximum value of the expression above.

For a given lambda the zero velocity curve describes the limit of where m is allowed. For a given lambda, an increasing rotation rate or decreased r₁, r₂ can be offset by increasing the amplitude of K.E. term - i.e. velocity rising. But where the rotation and G terms drop below the level of the constant (lambda term), then the kinetic energyda term), then the kinetic energy term would have to go negative - that is imaginary velocity, which is non-physical. So these areas are not allowed. Lambda will be defined by the initial overall energy of m.

We can find the zero velocity curves as follows:

Consider the motion of either body about the centre of gravity:

If we now take large values of lambda we see this can only be made (1) with very large x+y so that x²+y² is large (2nd and 3rd terms would be small). We then have an equation tending towards x²+y²=C i.e. a circle. It can also be satisfied with (2) x and y very small so the 2nd and 3rd terms are large (and the first is then very small).

The zero velocity curves are shown in the figure below as C subscripted with lambda and Ln where n is 1 to 3. The L3 case [case (1)] is the one with x and y large - i.e. the outside circles. Case (2) with x and y small corresponds to the L1 and L2 cases, circles centred around m₁ and m₂. That this is so we can see by looking at the 2nd and 3rd terms as x,y become small. The circle around m₁ is the case where:

so we tend to a circle:

centred at:

and similarly around m₂:

and the equation goes to:

centred at:

For smaller lambda the outer circle shrinks and the inner ones grow. Notice the "allowed" areas are outside the outer circle and inside the smaller ones. (Think what this means in terms of energy.) Eventually either the inner circles touch or one of the inner circles will touch the outer circle. In the case of the Earth-Moon system the inner circles touch, the limiting case being shown by the curve in the figure below labelled C with a subscript lambda and 1. In the inner curve the figure-of-eight has its crossover on the x axis. Note there is also a limiting outer curve as shown. With this value of lambda it is possible to go from the vicinity of m₁ to that of m₂ through the crossover (with zero velocity).

As lambda decreases further the crossover point L₂ turns into a neck. (Note this curve is not the trajectory of a vehicle or body - it merely defines the boundary to where m can go.)

Eventually as lambda gets smaller the case of curves C with subscript lambda and 2 is reached where inner and outer circles meet, initially at point L₃. By curves C with subscript lB>. By curves C with subscript lambda and 3 the "necks" are open at both ends and only the shaded areas are "not allowed" (energetically) to the body. These regions shrink as lambda gets smaller, and eventually end up as points L₄ and L₅. Fo lambda less than this limiting value nothing can be plotted, and the body m has so much energy it can move freely anywhere in the system. Small lambda means large initial energy. L₁, L₂, L₃, L₄ and L₅ are known as the Lagrangian points. The points L₄ and L₅ form equilateral triangles with m₁ and m₂.

... is to consider the balance of forces. These are all points of "equilibrium" where a body m would have zero velocity in this frame of reference. Thus at L₂ the centrifugal force on m plus the attraction of m₂ is "balanced" by the attraction of m₁. At L₁ the centrifugal force away from the centre of gravity is balanced by the attraction of m₁ and m₂. Similarly for L₃. These are points of unstable equilibrium in the m₁m₂ line, since we see any movement along the line destroys the balance of forces. (Example: take L₃. If the body m is moved to larthe body m is moved to larger x the centrifugal force rises while the attraction of m₁ and m₂ falls, so it continues to accelerate outwards. Movement towards m₂ means attraction to m₁ and m₂ rises and centrifugal force falls, so the body continues to accelerate towards m₂. Satisfy yourself that L₁ and L₂ are similarly points of unstable equilibrium.) Note that we have considered motion along the line m₁m₂ only. In the y direction there may be a certain degree of stability - and "halo" orbits have been used for spacecraft orbiting in the y-z plane about one of the L₁, L₂,L₃ points.

For the 3 Lagrangian points on the x axis we have:

and we can use this property to define an iterative methodefine an iterative method for solving for L₁, L₂ and L₃ by successive approximation.