Consider two gravitating bodies far from any other gravitating bodies or causes of other forces (but let's not assume either of them is negligible compared to the other). How do they move ?
Chose a frame of reference (for some neighbourhood of their entwined space-time trajectories) in which the sum of their momenta is entirely time-wards; i.e. their spatial momenta are equal and opposite.
Let the total mass of the bodies be M; let their several masses be a.M and
A.M respectively; a+A = 1. When they are a distance R apart, sub-divide the
line between them in the proportions a.R at A.M's end and A.R at a.M's end;
the point of sub-division is known as their
centre of gravity. Because
the total spatial momentum is zero, yet equal to M times the velocity of the
centre of gravity, the centre of gravity isn't moving; so take it as origin
for non-spinning spatial co-ordinates in conjunction with the time-like
co-ordinate choice which made the total (spatial) momentum zero. Then the
centre of gravity's trajectory in space-time is the time axis of the
co-ordinate system. I shall presume that M, A and a are constant; R, and the
direction of the line between the masses, will vary with time, t.
Let w be the result of
dividing a.M's velocity by A.R (i.e. a.M's distance from the origin); because
the centre of gravity (our spatial origin) always lies on the line from a.M to
A.M, dividing it in the same proportions, dividing A.M's velocity by a.R
(A.M's distance from the origin) must give −w, the same but negated
(because we measured both radii as positive, despite their being in opposite
directions; while the velocities encoded opposite direction by
fractional velocity, w, encodes both the fractional
rate of change of R, as its radial component, and the angular velocity of the
bodies about their common centre, as a circumferential component. Its
location on a diagram should be at
unit radius along the line from the
origin to a.M, extended past a.M if the
unit exceeds A.R; crucially,
it's on a.M's side of the origin; the corresponding vector to display on A.M's
side of the
unit circle would be −w.
Since a.M has velocity A.R.w, its momentum is a.A.M.R.w, while that of A.M will be −a.A.M.R.w. As a, A and M are constant, the rate of change of a.M's momentum is thus a.A.M.d(R.w)/dt; in accord with Newton's second law (for each action there is an equal and opposite reaction), that of A.M is the same but negated.
Newton tells us that the rate of change, a.A.M.d(R.w)/dt, of a.M's momentum is towards the origin with magnitude G.A.a.M.M/R/R, as long as the only force involved is their mutual gravitation (and neither R.w nor √(G.M/R) is near enough to light's speed to oblige us to worry about relativistic effects). The bodies' gravitational potential energies, each in the field of the other, are then −G.M.M/R times a.A.A and times A.a.a for a sum of −a.A.G.M.M/R.
Space-time has a metric, g; it's
constant and symmetric, so g(x,y)
= g(y,x), and I'll try not to be distracted by the fact that our co-ordinates
have been chosen to diagonalise it (with one component of one sign, three of
the opposite sign and all off-diagonals zero). Since w is the spatial
components of either body's velocity divided by the body's radius, it's
entirely spatial so g(w,w) only sees the spatial components of g, which I'll
presume to be positive-definite (making the time component negative).
Now, a.M's position vector divided by A.R is a unit vector, u. The velocity of a.M is A.R.w, yet is equally A.d(R.u)/dt = A.(u.dR/dt + R.du/dt), whence R.w = u.dR/dt + R.du/dt; so, in particular, w is in the plane spanned by u and du/dt. Now
and Newton tells us (see above) that this last is −u.G.M/R/R; i.e.
Introduce a unit vector, n, in the direction of du/dt, which is
orthogonal to u: this changes direction
in the same manner (subject to
the rotation which maps u to n) as u with dn/dt = −u.g(n,du/dt) and
du/dt = n.g(n,du/dt). We already know du/dt = w−u.g(u,w) and g(n,u) =
0, so we can infer g(n,du/dt) = g(n,w), whence
Our force-law is central, i.e. g(n,w.g(u,w)+dw/dt) = −g(n,u).G.M/R/R/R = 0 as g(n,u) = 0, so
making J = g(R.n,R.w) constant: multiplied by suitable masses, this gives the angular momenta of the bodies; a.M has angular momentum a.M.g(A.R.n,A.R.w) = a.A.A.M.J while A.M's is A.a.a.M.J and, as a+A = 1, the sum is S = a.A.M.J. As 1/(a.A.M) = (a.M+A.M)/(a.M)/(A.M) = 1/(a.M) +1/(A.M), J is the total angular momentum multiplied by the sum of inverses of the masses. [Thus S = J/(1/(a.M) +1/(A.M)), so 2.S/J = 2/(1/(a.M) +1/(A.M)), which is called the harmonic mean of the two masses.] J has dimensions area/time and, thanks to our choice of n's direction, J is non-negative.
In particular, g(n,w) = J/(R.R) can never change sign. Note that, for non-zero R, J = 0 iff g(n,w) = 0 iff w is parallel to u iff du/dt is zero iff u and n are constant. In any case,
Let Q be the dimensionless vector w.R.J/G/M −n, which is constant; we then have
so introduce the length L = J.J/G/M, the
(which is all Latin to me) of the
conic section (the curve of
intersection between a cone and a plane) our trajectory will turn out to be
following. Our trajectory can now be written as:
Discussion must, given that J is non-negative (by construction), now break into two cases: J = 0 and J > 0.
The following discussion involves the trigonometric functions sine and cosine; following habits from elsewhere in this web site, I'll use the Capitalised names Sin and Cos for the functions which map angles (in units of radians, turns, degrees or whatever) to scalars; and their lower-case equivalents, sin = (: Sin(x.radian) ←x :) and cos = (: Cos(x.radian) ←x :), for the closely related pure-scalar functions which satisfy cos' = −sin, sin' = cos and, in the complex plane, exp(i.z) = cos(z) +i.sin(z). The J=0 case will also exploit sinh(z) = (exp(z) −exp(−z))/2 and cosh(z) = (exp(z) +exp(−z))/2, making cosh(i.z) = cos(z) and sinh(i.z) = i.sin(z), with sin2 +cos2 = 1 = cosh2 −sinh2, sinh' = cosh and cosh' = sinh.
For the special case of circular orbit,
so, supposing the A.M body has radius C.R and density B, while the a.M body has radius c.R and density b, so that A.M = (4.π/3).B.C.C.C.R.R.R, a.M = (4.π/3).b.c.c.c.R.R.R, we find
enabling us to infer the densities – b from j, c and a, B from j, C and A. If we are observing the mutually-orbiting system from afar, possibly without knowledge of our distance from the system (hence without knowledge of the absolute value of R), we'll still be able to determine j (because 2.π/j will be the period of the mutual orbit, given ddu/dt/dt = −j.j.u), while the remaining quantities (A, a, C and c) are ratios between lengths, which we can infer from the angles we see, so insensitive to our ignorance of R. Thus we can determine the densities of mutually orbiting bodies from afar, at least when the trajectory is a circular orbit.
Furthermore, though one of the bodies may be too small to distinguish its radius from zero, if the other is large enough that we can see what fraction of their separation its radius is, then we can still determine its density, even though we cannot determine that of its partner. Thus Newton's contemporaries presumably had enough data to determine Jupiter's average density, if not especially accurately; at 1.33 g/cc it's only slightly denser than water, on average – less than a quarter of Earth's average density. This should at least have given some clue to the effect that Jupiter is a lot bigger than Earth. As moons of other planets showed up, the densities of Mars (3.94 g/cc, Earth / 1.4), Saturn (0.7 g/cc, Earth / 7.9), Uranus (1.3 g/cc, Earth / 4.2) and Neptune (1.76 g/cc, Earth / 3.1) must have come to light … oh, and Pluto (2.0 g/cc, Earth / 2.8), back when it was counted as a planet, but only relatively recently did we realize it was a binary system.
When J is zero, hence so is L and the above description of the trajectory becomes fatuous, 0 = 0, we still have Q = w.R.J/G/M −n = −n constant, so conclude that n and u are constant: the bodies are confined to a fixed line through the origin, with no angular velocity; their momenta are collinear and there is no mutual angular momentum. Now, with u and n constant, w is simply u.g(u,w) = u.dR/dt/R. Our dynamical equation reduces to:
Multiplying this last by a.A.M/2 yields the total energy; dividing it by 2.G.M yields a constant with dimensions 1/length; call this K; it might be zero, positive or negative. We obtain ±dR/dt = √(2.G.M.(K+1/R)) so integrating dt = ±dR/√(2.G.M.(K+1/R)) will complete our solution.
At this point we can examine energies: the kinetic energy of a.M is (1/2).a.M.(d(A.R)/dt)2 = (1/2).a.A.A.M.(dR/dt)2 and that of A.M is (1/2).A.a.a.M.(dR/dt)2 for a sum of (1/2).a.A.M.(dR/dt)2 = (1/2).a.A.M.2.G.M.(K+1/R); as observed above, the total gravitational potential energy is −a.A.G.M.M/R; adding this to the kinetic energy we get a.A.M.((dR/dt)2 −2.G.M/R)/2 = G.(a.M).(A.M).K which is constant; it's the gravitational potential energy of either body in the field of the other at separation 1/K. So the (constant) total energy for the collision course is:
Now, to solve (dR/dt)2 = 2.G.M.(K+1/R), consider three branches, depending on the sign of K:
with K+1/R > 0, hence R < −1/K, so introduce an angular parameter s with −K.R = (1 −cos(s))/2 = (sin(s/2))2 so that
then, choosing t = 0 at a moment when s = 0, hence R = 0, the given constant is zero so we have sin(s) −s = ±2.K.t.√(−2.G.M.K) with R = (cos(s) −1)/2/K describing our solution.
Now, R = 0 happens when cos(s) is 1, i.e. when s.radian is a whole number of turns. Since sin(s) −s changes by 2.π between such moments, at each of which sin(s) is 0, t changes by −π/K/√(−2.K.M.G), giving this as the time it takes the bodies to fly apart after one collision, then return to collide once more. Of course, in practice, the bodies shall interact in some non-gravitational manner when they collide, so you'll only see part of one cycle of this process; the analysis ignores this, effectively pretending that they either pass through one another, unaffected, or bounce off one another perfectly elastically.
Writing H for −1/K, the maximum value R attains, we get a total energy of a.A.G.M.M/H and sin(s) −s = ±2.(t/H).√(2.G.M/H) with a period of H.π.√(H/2/G/M).
By taking t = 0 when R = 0, at the moment of collision, this becomes
for which the bodies fall together, pass through or bounce off one
another (the given solution doesn't tell you which, or care) at t=0 then fly
apart in the time-reverse (possibly
centrally inverted spatially) of
their earlier collision trajectories. Both long before and long after the
collision, the bodies are very far apart and essentially stationary; in the
limit as long tends to infinitely, they are infinitely far apart and actually
By analogy with the negative case, consider a parameter s for which K.R = sinh(s/2)2 = (cosh(s) −1)/2 so that
again – choosing t = 0 when s = 0, hence R = 0 – we can make the constant zero, giving sinh(s) −s = ±2.K.t.√(2.G.M.K) while R = (cosh(s) −1)/2/K. For large values of s, cosh(s) and sinh(s) are roughly equal and s (despite being large) is negligible compared to them, yielding R = t.√(2.G.M.K). This solution describes bodies intially falling towards one another at relative speed √(2.G.M.K) and ultimately flying apart at the same relative speed.
Writing V = √(2.G.M.K) for this initial and final relative speed, yielding 2.K = V.V/(G.M), we find that the solution is sinh(s) −s = ±2.V.V.V.t/G/M.
The similarities between the K<0 and K>0 cases arise from the broken duality between [sin, cos] and [sinh, cosh]. Still, enough of the unrealistic special case with exactly no angular momentum …
If J = √(L.G.M) is non-zero (in which case our choice of n makes it positive), we can divide G.M.Q = J.w.R −G.M.n by J to obtain a constant velocity V = G.M.Q/J = w.R −n.G.M/J = w.R −n.J/L, giving Q = J.V/G/M = V.L/J, so expressing our trajectory as
If ever u is parallel to V, g(n,V) will be zero, making R = L; i.e. L is the radius at which the trajectory crosses the line, though the origin, in the direction of V. [V is perpendicular to the focal axis of the conic that R.u traces out; V is of order 36 km/s for Earth's orbit around the Sun; one thousandth of G.M/J]. In any case, we now have
and R.w, our bodies'
normalised velocity, consists of a
constant, V, plus a rotating
cirumferential component of constant
magnitude, J/L. Any resemblance to the
epicycles of the ancient Greeks
can be construed as evidence they were less wrong that they're sometimes
accused of ;^)
In the following analysis, I'll make use of a constant speed, v, which is
V's magnitude, i.e. v is a positive scalar with v.v = g(V,V). Likewise, let e
be the magnitude of our dimensionless vector Q = V.L/J, i.e. e = v.L/J,
yielding v.v = e.e.G.M/L. This e is known as the
eccentricity of the
trajectory. Since n was so chosen as to make J non-negative, and we're dealing
with the case of non-zero J, e cannot be negative.
The kinetic energy of a.M is a.M.g(A.w.R,A.w.R)/2, that of A.M is A.M.g(a.w.R,a.w.R)/2, for a sum of a.A.M.g(w.R,w.R)/2, so
in which only the last term varies; and its value, a.A.G.M.M/R, is just −1 times the total potential energy, so the sum of potential and kinetic energy is
Thus the total energy is positive iff e.e > 1, i.e. the eccentricity is bigger than 1. When R = L, the kinetic energy is a.A.(e.e+1).G.M.M/L/2, which is −(e.e+1)/2 times the potential energy.
The rotating frame of reference induced by n and u sees w.R = V +n.J/L as
constant circumferential velocity n.J/L plus a
offset, V, which adds a sinusoidally varying component g(n,V) to the angular
velocity and contributes a complementary sinusoidally varying component g(u,V)
as the radial velocity; but note that the sinusoid's angular velocity is the
angular velocity to which g(V,n) is contributing, so the sinusoid's phase
doesn't vary linearly with time.
Introduce an angle θ specified by: g(n,V) = v.Cos(θ), g(u,V) = v.Sin(θ) – i.e. chose (positive) y-axis parallel to V and look at the plane of the trajectory from whichever (spatial) side makes the angular velocity positive (anti-clockwise), to infer the x-axis (a quarter turn clockwise from the y-axis), then use the usual polar co-ordinate measured anti-clockwise with θ zero on the positive x-axis.
Since V = n.v.Cos(θ) +u.v.Sin(θ) is constant, so are both its magnitude, v, and its direction, V/v = n.Cos(θ) +u.Sin(θ). Our trajectory is now expressed as
with dR/dt = g(u,V) = v.Sin(θ). We can now re-write the total kinetic energy as
If we write X = R.Cos(θ), Y = R.Sin(θ), we now have L = R + e.X whence
which describes an ellipse for e.e<1 and a hyperbola for e.e>1. In either of these cases, substituting X=0 yields Y.Y = L.L so Y = ±L (as already anticipated); while substituting Y=0 yields X = L.(±1 −e)/(1−e.e) = L/(e±1). For the case e.e = 1, we get Y.Y = L.L −2.L.X, or 2.L.X = L.L − Y.Y, which is a parabola.
Our dynamical equation now reduces to θ's time-dependence,
so time is the integral of L.L/J/(1+e.Cos(θ))2 with respect to θ/radian, plus some arbitrary offset – albeit this integral is no picnic ! Again, we break into three cases, according as the total energy is positive, zero or negative; which is to say, according as e is greater than, equal to or less than 1. The X,Y analysis above revealed that these cases give hyperbolic, parabolic and elliptical orbits, respectively.
From L/R = 1 + e.Cos(θ) with L, R and e positive, we can infer Cos(θ) > −1/e. This gives us a wedge-shaped portion of the plane from which the (R,θ) representation of our bodies' trajectories are excluded; the apex of the wedge is our origin and its two bounding rays point outwards in the directions with Cos(θ) = −1/e; these are in the −ve X half-plane, one in its +ve Y quadrant, the other in its −ve Y quadrant. The portion between these, straddling the −ve X-axis, is the wedge-shaped region from which the orbit is excluded.
Our hyperbola is
centred on the point X = e.L/(e.e−1), Y = 0,
i.e. R = e.L/(e.e−1), θ = 0. If we translate the exclusion-wedge
in the X-direction through e.L/(e.e−1), the hyperbola lies entirely
within the translated wedge, though outside the un-translated wedge; it is
confined to the >-shaped trench between the two wedges' boundaries. In its
(R,θ) representation, our bodies' trajectories come from far out along
the −ve Y branch of the translated wedge's edge; move steadily towards
the apex of the > and away from this leading edge, speeding up as they move
inwards; cross the line θ = −turn/4 at R=L; swing round through
the angle in the trench, crossing θ=zero at R=L/(1+e) < L/2; then
escape in like manner, though time-reversed and reflected in the X-axis, to
ultimately approach the translated wedge's edge's +ve Y branch.
The relative speed of the two particles while very far apart is just v.√(1−1/e/e) = v.√(e.e−1)/e = J.√(e.e−1)/L. The maximum relative speed is J.(e+1)/L, which happens when θ=zero, R=L/(1+e). The ratio of these two speeds, dividing the former by the latter, is √((e−1)/(e+1)); it tends to 1 for large e but tends to zero as e approaches 1. When far from each other, the bodies' velocities are almost exactly either towards or away from one another; the magnitude increases with e; near e = 1, when the internal angle of the wedges is tiny (and the outer vertex is far from the origin, to keep the Y-wards width still L), they are barely moving (relative to one another).
This will swing round the
origin in very similar manner to the hyperbola, from far away in the −X
direction, through the −ve X, −ve Y quadrant, crossing into the
+ve X, −ve Y quadrant at R = L on the −ve Y-axis, swinging through
this quadrant to pass through the X-axis at X = L/2 and a top speed of 2.J/L
before followinng a time-reversed mirror image of its infalling
trajectory. Its speed tends to zero at large radius. This can be thought of
as the limiting case of the hyperbola as e tends to 1 from above, albeit the
wedge has become the whole plane (its outer apex is infinitely far away and
the only place excluded is the negative X-axis); θ's value is
arbitrarily close above −turn/2 for large negative t and tends to
+turn/2 from below as t tends to infinity. It can also be thought of as the
limiting case of an ellipse (see below) with the distant end of the major
We can exploit 1+Cos(θ) = 2.Cos(θ/2)2 to obtain J.dt/L/L = 4.dθ/Cos(θ/2)4/radian. Substituting z = Tan(θ/2), we obtain 2.radian.dz/dθ = Sec(θ/2)2 = 1+z.z whence
whence z.(1 +z.z/3) = J.t/L/L/8 = 3.t.(G.M)2/(2.J)3 plus a constant which we can make zero by taking t = 0 when z = 0, i.e. when θ is zero. We start with z large and negative (θ just increasing from −turn/2) and end with z large and positive (θ tending towards turn/2). For large t, Tan(θ/2) = z is then approximately half the cube root of 3.J.t/L/L or, equally, approximately 1/2/J times the cube root of 3.t.(G.M)2.
with (1−e.e).Y.Y +((1−e.e).X +e.L)2 = L.L, R.u's orbit encloses an area π.L.L/(1−e.e)/√(1−e.e). The rate at which the R.u ray sweeps out area is g(R.n,R.w)/2 = J/2, which is constant, and this ray sweeps out the entire area of the ellipse in each orbit, so the period of the orbit is
Substituting J = √(L.G.M) into this gives a formula in which L/(1−e.e) appears with power 3/2; and L/(1−e.e) is just the semi-major axis of our ellipse, since L/(1−e) +L/(1+e) = 2.L/(1−e.e). Thus, writing H for the semi-major axis, we get period 2.π.H.√(H/G/M) regardless of eccentricity. Furthermore, E = (e.e−1).G.a.M.A.M/L/2 = −a.A.G.M.M/H/2, so (for a given pair of bodies) both the total energy and the period depend only on the semi-major axis of the orbit. The product of energy and period, −π.a.A.M.√(G.M.H) has the same dimensions as angular momentum, S; so introduce
The orbit is circular when e is 0, making f = π so −E.period = π.S with L equal to the semi-major axis (which, in this case, is the radius).
Since R, u and n are expressed in terms of θ, relating θ to time (i.e. integrating up the equation dt.J/L/L = dθ/radian/(1+e.Cos(θ))2) is all that remains for a full, general solution of Newton's two-body problem.
Substitute u = k.Tan(θ/2), for some constant k to be chosen in due course: this gives Cos(θ) = (k.k −u.u)/(k.k +u.u), du = k.Sec2(θ/2).dθ/2/radian, so dθ/radian = 2.du.Cos2(θ)/k and we get
whence, picking t = 0 at some moment when θ is 0 or a half turn, hence u is zero,
The final du is easy; in the other two terms, substitute k.Tan(v).√(1+e) = u.√(1−e), i.e. v = arcTan(Tan(θ/2).√((1−e)/(1+e))), so that du = k.Sec2(v).(dv/radian).√((1+e)/(1−e)) while k.k.(1+e) +(1−e).u.u = k.k.(1+e).(1 +Tan2(v)) = k.k.(1+e).Sec2(v), to get:
in which Cos2(v).dv/radian = (1 −Sin2(v)).dv/radian = dv/radian +d(Cos(v)).Sin(v) = dv/radian +d(Cos(v).Sin(v)) −Cos(v).d(Sin(v)) = dv/radian +d(Cos(v).Sin(v)) −Cos2(v).v/radian, the last term in which is just what we started with negated, hence 2.Cos2(v).dv/radian = dv/radian +d(Cos(v).Sin(v)), so:
which should at least let us compute t from θ, albeit going the other way is apt to be rather tricky !
Note that, since the energy (for a given pair of bodies) and the period only depend on the semi-major axis of the orbit, in each case as a monotonic function of H, we can infer that orbits which pass through a given point at equal speeds must necessarily have equal periods – and vice versa. (I first met this detail as a school-boy in the form of a question specifying that an object, instantaneously at rest in a gravitational field, is blown up by a bomb which causes every part of the object to get a share of the released energy, in proportion to its mass; which amounts to specifying equal speeds: all the fragments shall return to their start-point at a common later moment – ignoring any which run into something else along the way.)
The period of anything is an amount of time per cycle; so, properly, the
period of our elliptical orbit (J>0, e<1) should be construed as an
rate of orbit of
2.π.L.L/J/(1−e.e)/√(1−e.e) per turn, hence
equally as L.L/J/(1−e.e)/√(1−e.e) per radian. This turns f
into E/S/(rate of orbit) and makes 1/f equal to 2.(1−e.e).radian unless
we give S units that include a 1/angle term to cancel the angle term arising
here in f. The proper rôle and place of units of angle, in particular
in connection with
angular momentum, is an area that causes me some
doubts: notably, Planck's constant is an energy times period, like f.S, so
should really involve the unit of angle – with the result that
multiplying it by turn yields the quantity usually known as Planck's constant,
while multiplying it by radian yields Dirac's constant, which is a factor of
Note, also, that L.L/J/(1−e.e)/√(1−e.e)/radian = J.J.J/G/G/M/M/(1−e.e)/√(1−e.e) =
which encourages a review of the above with attention to J/√(1 −e.e) or L/( 1−e.e) as the terms to which to pay attention.
In all of the above, I have used a non-rotating frame of reference to
describe the system; however, since the elliptical orbit has an over-all
average angular velocity, it makes sense to
subtract off this rotation
and describe the orbit from the point of view of
a rotating frame of reference. The
orbit (or inverse period) alluded to above is
J.radian.(1−e.e).√(1−e.e)/L/L, so this is the angular
velocity to use for our rotating frame. The resulting frame will be useful
for considering the motions of smaller bodies in the
neighbourhood of such a mutually orbiting pair.
Our natural point of departure from the above is the pair of equations
and we are given that J>0 and 0≤e<1. We are changing to a co-ordinate system in which R is unchanged but θ is replaced by
so that, introducing D = L.L/(1−e.e)/√(1−e.e) for the result of dividing our ellipse's area by π,
For circular orbit, e = 0, we have R = L and φ constant
I think the right approach will be to assume small e, derive the dynamical equations in the rotating frame, apply some suitable approximations and solve by perturbation – expect another epicycle ;^)Written by Eddy.