The Stability of the Double Binary System PH1


In October of 2012, the first four-star planet was confirmed by the Planet Hunters program from Yale. The planet – called Planet Hunters 1 (PH1) – is a circumbinary planet meaning it is orbiting a pair of stars instead of just one. Furthermore, orbiting that pair of stars is another pair meaning this planet is in a system with four stars total. With so many large forces acting on this system, the stability is, of course, a question of concern. According to the paper announcing PH1’s discovery, though, “the system is indeed stable on gigayear timescales” (Schwamb 2013). This implies there are formation possibilities never considered and is inspiring further study of PH1.


Using the Kepler Space Telescope, astronomers and planet hunters have discovered and confirmed 1,020 planets from the time of its launch in 2009 and the time of this writing. There are still over 4,600 extraterrestrial objects that seem to be planets; they just have not been studied enough to be confirmed yet. Along with those confirmed and unconfirmed planets, Kepler has discovered 2,165 binary star systems({Kepler Website}). The question then becomes, though, are these binary systems stable? The goal of this paper is to

  1. Analyze Holman and Wiegert’s paper Long-Term Stability of Planets in Binary Systems (Holman 1999) to fully understand what makes such a system stable and

  2. Apply theses finding to the system PH1 to find if it, too, is stable

The Simulation Process

Holman and Wiegert begin by addressing the question of where planets can last for long periods of time instead of asking how they are formed. First, they must consider what size planets can persist for long periods of time. As Szebehely found, as the mass parameter, \[\mu= \frac{m_2}{(m_1 + m_2)},\] increases from \( \mu=10^{-6} \to \mu \approx 0.3\) the dimensionless radius for stability of the orbiting planet increases exponentially from \(r=1 \to r=2.4\) where the radius quickly decreases (Szebehely 1984). This means the more similar the two stars are in mass, the larger the planet must be to remain stable. Knowing that, Holman and Wiegert then describe their simulation. They use an elliptic restricted three body system where the planets are test particles that do not interact with each other so the simulation can be done with multiple varying planets at the same time. During simulation, they look for “close encounters” — which they define as a passage within 0.25 of binary separation — and for escape orbits. So any test particle that comes too close or too far away from the stars are deemed unstable and removed from the integration. They then determine what they call the critical semimajor axis which is the “semimajor axis at which the test particles at all initial longitudes survived the full integration time.” They ran this simulation for many orbit shapes varying from perfectly circular to eccentricities of 0.8.

Simulation 1

For the first simulation, Holman and Wiegert used \(\mu=0.5\) and \(e=0.5\) for the inner region which means the planet was orbiting only one of the stars. After this simulation they found a critical semimajor axis of \(a_c = 0.12\). Their results showed that the particle was more stable the smaller the initial semimajor axis was. At a starting semimajor axis of \(a=0.17\), no planet survived more than 30 binary periods. Then, as the semimajor axis decreased, the planets lasted more binary periods. The planet closest to its parent star lasted longer than any of the other stars with the same semimajor axis for every trial but one. Interestingly, though, the increase in binary periods lasted is very slow until a semimajor axis of \(a=0.14\) where it jumps to a median binary periods lasted of 175.5. Then, at the critical semimajor axis, there is a huge jump where every planet lasts 10,000 binary periods. This implies a couple of things; first, the distance from the star is not as influential in the stability as the semimajor axis. Second, the critical semimajor axis is a threshold. The semimajor axis does not gradually become more stable, it is either stable or not.

Comparison with Earlier Work

The authors then compare their results to that of Rabl and Dvorak who did an almost identical simulation 11 years prior (Rabl 1988). One major difference was simply the computing power. Rabl and Dvorak only ran the integration for 300 binary periods because computers were slow and expensive meaning running the simulation longer would have been uneconomical. Even with the lower precision, however, this simulation found very similar results. Looking at Rabl and Dvorak’s congruent table (\(e=0.5\)), the same patterns arise: as the semimajor axis decreases, the stability increases. The planet closest to the stars is generally more stable compared to further planets with the same semimajor axis. And, most importantly, at \(a=0.12\), every planet remained stable for more than 300 binary cycles. Holman and Wiegert calculated a best fit line for their critical semimajor axes which came to \[a_c = .274-.338e+.051e^2 ,\] which they compared to Rabl and Dvorak’s: \[a_c = .262-.254e-.06e^2 .\]