One of the most beautiful results in dynamics is the thermal eccentricity distribution. Suppose we have a population of binary stars which is, in some sense, “thermalized.” That is, these binaries have all interacted with each other and exchanged energy many times and have reached statistical equilibrium. We would like to know what these binaries look like. In particular we would like to know the distribution of their periods and eccentricities. The derivation of this result is due to the British astronomer J. H. Jeans in a 1919 paper and this post loosely follows Jeans’s derivation. I have updated some of the notation, simplified the derivation somewhat, and provided some motivation behind Jeans’s sometimes cryptic coordinate transformations.

It is clear from the outset that the problem as we have posed it has a small issue. It is not obvious that there is any statistical equilibrium that can be reached at all. When these binaries interact, it seems that some of them will dissociate and others of them will form triples. So a thermalized population of binaries won’t be a population of binaries at all, but a population of binaries mixed with single stars and triples (and given that these objects will all interact with each other, there will also be higher order systems as well). If this is the case, then there can be no such thing as a thermal population of binaries. This objection turns out to be valid because a population of binaries will naturally form single and triple systems. But we will here make a small swindle and suppose that such a thermalized population of binaries exists. This population need not have to have come about naturally, but we can imagine that it was instead created ab initio and then ask what its properties will be knowing that there is equipartition of energy. We should not object to this swindle too much because Jeans is already notorious for a far greater swindle. Moreover, we will find that, having made this swindle, the distribution of eccentricities takes a remarkably simple form and the distribution of periods validates the very objection we had raised!

So let’s begin. By assumption, we have a population of binaries which is thermalized. This is to say that the distribution of energies follows a Boltzmann distribution:

where \(T\) is some measure of the thermal content of the system. Now, the energy of an individual binary system is

where \(r\) and \(v\) are the relative distance and velocity, respective, \(M\) is the total mass of the system, and \(\mu\) is the reduced mass of the system,

The total number of systems in a differential element of phase space, \(d \tilde{V}\), is

where \(d\Omega\) is a differential solid angle. It helps if we decompose the velocity into its radial and tangential components: \(v^2 = \dot{r}^2 + \tau^2\). We can pick \(\dot{r}\) to be along the \(z\) direction, in which case we have \(dx \, dy = \tau \, d\tau \, d\phi\), where \(\phi\) is an azimuthal angle. Now, integrating over \(\Omega\) and \(\phi\), which gives us a factor of \(4\pi\) and \(2\pi\), respectively, we get

We so far have been using real, physical coordinates to describe the binaries. But what we would like is to transform this distribution function into a set of coordinates which represent the overall properties of the binary, like its specific energy and angular momentum. As it turns out, using the angular momentum itself complicates the derivation quite a bit and the inverse of the angular momentum is the better coordinate to use. So we have two new coordinates:

The shape of an orbit is fully determined by its energy and angular momentum (recall that we have already integrated over all possible orientations). We are therefore missing only one coordinate, namely a phase angle specifying where the stars are in their orbits. For this purpose we use the true anomaly, \(\theta\), which is related to \(r\) by

where \(A\) is the magnitude of the Runge-Lenz vector:

We may rewrite the equation for the true anomaly in this simple form:

To make this transformation we need to calculate the determinant of the Jacobian. Thanks to Jeans’s judicious choice of coordinates, this turns out to be simple by noting that neither \(k\) nor \(\theta\) depend on \(\dot{r}\), so

This implies that the volume element transforms as

Writing the old, physical coordinates in terms of the new ones, we have

and

Now recall the relationship between the magnitude of \(A\) and the energy and angular momentum:

This means that we can rewrite \(\dot{r}\) as

Putting all this into the distribution function, we find

We can now integrate this over \(\theta\) to get the distribution function in terms of the variables we wanted all along—the energy and angular momentum.

After performing the integral we find the distribution function to be

Here things start to simplify dramatically. First, note that \(A\) is related to the eccentricity by

Substituting into the distribution to rid ourselves of \(A\) forever, we get

Of course, we’re not so much interested in the distribution of angular momenta as we are in the distribution of eccentricity. So now that the distribution is fairly simple, let’s change coordinates once more, this time moving from \(k\) to \(e\). The two are related by

When we make this coordinate transformation we find that the distribution function has simlpified dramatically to

The eccentricity dependence of the distribution function is contained entirely in the term \(e \, de\) and is independent of the energy. We may therefore write the distribution function in terms of the eccentricity alone as

where we have now properly normalized \(f\). After all our work, this is our grand result—the thermal eccentricity distribution. If a population of binaries has thermalized, the eccentricities will be distributed linearly and nearly circular orbits will be rare. In fact, the median eccentricity of this distribution is \(1/\sqrt{2}\).

Let us turn now to the energies of the binaries. Here the result seems elegant:

However, when we try to normalize this distribution function we encounter a problem—it diverges as \(E \to 0\) and as \(E \to -\infty\). This result would appear to contradict our entire analysis—a thermal population of binaries cannot exist at all! And indeed, this is the case. The divergence of this distribution function at both limits is a direct consequence of Heggie’s law. Heggie’s law, first stated in the dense monograph on three-body dynamics by the British astronomer Douglas Heggie, states that in general, when a passing star interacts with a soft binary (i.e., a binary whose binding energy is less than the kinetic energy of the incoming star), the binary will tend to become softer and when a passing star interacts with a hard binary (i.e., a binary whose binding energy is greater than the kinetic energy of the incoming star), the binary will tend to become harder. This means that if we start with a population of binaries of all different sizes and moving with a variety of velocities, after many interactions, the hard binaries will have hardened indefinitely (they will have \(E \to -\infty\)) and the soft binaries will have softened indefinitely (i.e., disrupted). So our inability to normalize the distribution function actually fits in quite well with our understanding of dynamics.

Unfortunately Jeans predated the discovery of Heggie’s law by many decades, so he was unable to interpret the energy distribution as easily as us. Jeans nevertheless noted three conclusions from his derivation:

1. “There would be no correlation between period and eccentricity.”

2. “The eccentricities would be distributed according tot he law \(2e \, de\)—in other words, all values of \(e^2\) would be equally likely.”

3. “The periods would conform to the law of distribution.”

Jeans then compared these results with the observations available to him in 1919 to see how thermalized the population of binaries that had been observed was. It is obvious from the distribution of energies that the binaries cannot be thermalized. But it is possible (and indeed it is the case) that the eccentricities could thermalize before the energies. Nevertheless, the sample of binaries Jeans uses do not follow a thermal eccentricity distribution—there is a deficit binaries with eccentricities above 0.6. Moreover, the sample of binaries also exhibits a strong correlation between period and eccentricity. Jeans therefore concludes that binaries have not yet thermalized, and so the distribution of their orbital parameters must provide information about their formation.

Jeans’s sample (which he draws from The Binary Stars by Robert Aitken) was not ideal. It’s small (87 stars) and biased. How does a modern sample hold up? An excellent overview of the properties can be found in the review by Duchene & Kraus (2013). In short period binaries the eccentricities are almost all very low due to tidal circularization. In longer period binaries (\(P \gtrsim\) 100 days) the eccentricity distribution is consistent with flat, although there appears to be a dearth of systems at very high eccentricities. At any rate, modern samples confirm Jeans’s original observation that binaries in the field do not appear to be thermalized. Nevertheless, the modern samples do not find a correlation between eccentricity and period aside from the circularization of very short period binaries already mentioned.

It is not very surprising that binaries in the field do not exhibit a thermal eccentricity distribution because even over the lifetime of the Galaxy nearly all binaries have not had time to interact with each other even once. In a more compact system like a globular cluster we would expect the system to be more thermalized because the binaries would have many opportunities to interact with one another. Unfortunately it is generally difficult to measure the orbital properties of binaries in globular clusters—even getting an accurate measurement of the binary fraction has proven to be difficult. It is possible to get accurate orbital periods from millisecond pulsars and these seem to show a bias towards circular orbits. But it is reasonable to believe that millisecond pulsars represent a special case because the partner would necessarily have to be very close to the neutron star to spin it up and therefore a bias towards circular orbits would be expected. Without observations, we have to resort to numerical simulations, which seem to indicate that binaries in globular clusters should exhibit a thermal eccentricity distribution, but this cannot be concluded too definitely because it is standard for numerical simulations to begin with a thermal eccentricity distribution. Nevertheless it is a good confirmation of the theory that the thermal eccentricity distribution is maintained.

The thermal eccentricity distribution was also derived by the Soviet-Armenian astronomer Victor Ambartsumian in 1937 and we have a separate post about his derivation here. Ambartsumian’s derivation is more general and shows that the same eccentricity distribution can occur even if the binaries are not thermalized. For this reason the thermal eccentricity distribution is sometimes referred to as the Ambartsumian distribution instead.