The Dynamics and Stability of Circumbinary Orbits
Abstract
We numerically investigate the dynamics of orbits in 3D circumbinary phasespace as a function of binary eccentricity and mass fraction. We find that inclined circumbinary orbits in the ellipticallyrestricted threebody problem display a nodal libration mechanism in the longitude of the ascending node and in the inclination to the plane of the binary. We (i) analyse and quantify the behaviour of these orbits with reference to analytical work performed by Farago & Laskar (2010) and (ii) investigate the stability of these orbits over time. This work is the first dynamically aware analysis of the stability of circumbinary orbits across both binary mass fraction and binary eccentricity. This work also has implications for exoplanetary astronomy in the existence and determination of stable orbits around binary systems.
keywords:
celestial mechanics — stars: binaries — planetary systems1 Introduction
The recent discovery of a circumbinary disk around the microquasar SS433 (Blundell et al., 2008) along with evidence that such disks may be dynamically coupled to the accretion and outflow from such systems (Artymowicz & Lubow, 1996; Regös et al., 2005; Doolin & Blundell, 2009; Perez M. & Blundell, 2010) has led us to investigate the behaviour of orbits encompassing binary systems.
Additionally, in this era of exoplanetary astronomy we are finding new planets almost every week, but only a minority those so far discovered are in binary systems (e.g. Lee et al., 2009; Beuermann et al., 2011; Qian et al., 2011). With an unknown but significant fraction of stars confined to binaries, and with new technologies and methods to detect circumbinary planets (e.g. Schwarz et al., 2011), it is crucial that we understand the dynamics and stability of circumbinary orbits.
In the following sections we investigate a circumbinary nodal libration with the aid of finely timesampled numerical studies. We compare our results to analytic work performed by Farago & Laskar (2010) and investigate the accuracy and limits of their model.
With a second suite of numerical simulations we then determine the stability of the librating circumbinary orbits as a function of binary eccentricity and mass fraction.
2 Methods and terminology
2.1 Orbital Elements
A general orbit in 3D is described intuitively by the Keplerian orbital elements (, , , , , ) illustrated by Figure 1. The eccentricity , semimajor axis and true anomaly describe the motion of a body in its orbital plane, whilst the inclination , longitude of the ascending node and argument of perihelion describe the orientation of the orbital plane with respect to some reference plane and direction. The relationship between these orbital elements and a body’s state vector (position and velocity) may be found in Green (1985) amongst others.
In this coordinate space all trajectories are represented uniquely, with a Keplerian orbit having the property of conserving all quantities but the true anomaly , which describes the exact position of a body on its orbital path.
In the circumbinary case we define the reference plane as the natural plane of the internal binary, and the reference direction as the vector along the binary’s line of apses
2.2 Numerical setup
We have performed 3D numerical simulations of initially circular () circumbinary
We apply a customised adaptive stepsize fourth and fifthorder RungeKutta integrator to integrate suites of test particles as a function of binary eccentricity and mass fraction where .
We track the orbital elements of each test particle about the centre of mass of the binary during integration, outputting timelapsed snapshots to a database. We take advantage of the speed, organisation and SQL functionality of the database to study and accurately fit to the behaviours of orbits which we present in further sections.
Each test particle is also monitored for instability. Unstable orbits are identified and removed during integration where a test particle is perturbed sufficiently from its initial orbit to approach either star, or if it evolves onto an unbound trajectory (). This paper is not concerned with the ultimate fate of any test particle that experiences a close encounter with either stellar body. These test particles are rejected from the simulations and the precise scattering is not computed.
Postsimulation stability criteria are applied to identify and reject test particles which do not quite reach escape velocity.
3 The nodal libration
3.1 Suite of simulations
Our first suite of simulations was designed to be extensively timesampled to expose the dynamics of circumbinary orbits. The binary eccentricity and massfraction parameter space that we explore is laid out in Table 1 whilst the initial phasespace sampling of test particles is laid out in Table 2. One additional spherical shell of test particles at a semimajor axis was also integrated for .
Orbital Element  min  max  

eccentricity  0  0.6  0.1 
mass fraction  0.1  0.5  0.1 
Orbital Element  min  max  

semimajor axis  1.5  10  0.5 
inclination  
longitude of the ascending node  
true anomaly  
simulation length and snapshot 
= binary semimajor axis and = binary orbital period.
3.2 Description
As a typical test particle is integrated over the course of a simulation its semimajor axis and eccentricity remain constant. We discover a nodal libration mechanism in the inclination and longitude of the ascending node of all inclined orbits.
This is best visualised in a polar slice through parameter space, which may be termed a surface of section, as illustrated by Figure 2. Each line in Figure 2 is traced out by an individual test particle over the course of a simulation. Figure 2 reveals four distinct populations that we now examine in turn.
Prograde orbits (green)
An inclination corresponds to a coplanar prograde circumbinary orbit. Orbits of inclination , whilst not necessarily coplanar, we refer to as prograde. The prograde region of phase space therefore extends out from the centre () of the surface of section (Figure 2) and encompasses all orbits of a similar behaviour.
Prograde orbits exhibit a precession in the longitude of the ascending node . This evolution in produces clockwise paths around the surface of section shown in Figure 2.
Retrograde orbits (blue)
An inclination corresponds to a coplanar retrograde circumbinary orbit. Orbits of inclination , whilst not necessarily coplanar, we refer to as retrograde. The retrograde region of phase space therefore extends inwards from the outer limit () of the surface of section (Figure 2) and encompasses all orbits of a similar behaviour.
Retrograde orbits exhibit a precession in the longitude of the ascending node . But counter to the prograde orbits the retrograde evolution in produces anticlockwise paths around the surface of section (Figure 2).
We expect that the precession in observed in closetocoplanar prograde and retrograde orbits is due to a coupling between the specific angular momentum of test particles on inclined orbits and the angular momentum of the binary. Such a coupling would exert a torque on the test particle producing a precession in the ascending node, akin to gyroscopic precession.
Islands of libration (red & purple)
An inclination corresponds to a circumbinary orbit which is exactly perpendicular to the binary plane. Figure 2 shows two very clear libration islands centred on , . A test particle on an orbit within a region of libration has its inclination and ascending node coupled to precess about the centre of libration. For both regions of libration this precession is anticlockwise.
3.3 The geometry of the surface of section
The geometry of the surfaces of section shown in Figure 2 reveal a dependence on the internal binary eccentricity. More specifically, the extent of the two regions of libration can be seen to scale with binary eccentricity. A circular binary exhibits no libration islands, whereas for a binary of the libration mechanism is becoming the dominant behaviour in phase space. We quantify the extent of the libration regions in § 3.5.
Inspection of surfaces of section across values of binary mass fraction and radius lead us to conclude that the geometry is both mass fraction and radius independent. The period of the precession of each test particle however does show a strong dependance on and , which we explore in § 3.7.
Kozai cycles
The geometry of the polar (,) surface of section (Figure 2) with its islands of libration appears similar to that of the Kozai mechanism. Kozai (1962) showed analytically that an inclined circumstellar
Symmetry
The surface of section (Figure 2) is a 2D projection of the surface of a unit sphere, where corresponds to the polar angle and to the azimuthal angle. The central point of the surface of section () corresponds to a coplanar prograde orbit (see Figure 1) and at this point the longitude of the ascending node is undefined. Equivalently one may regard this point as situated at the north pole of a sphere, with the south pole at (a coplanar, retrograde orbit). Points of identify the equator. This type of spherical projection is known in cartographic circles as the azimuthal equidistant projection.
This unit sphere essentially defines the direction of the specific angular momentum vector of a test particle. We note three particular planes of symmetry through this sphere. These may be specified by considering components of a circumbinary orbit’s specific angular momentum along axes (i) parallel to the line of apses of the binary, (ii) perpendicular to the line of apses of the binary
green:
blue:
red:
3.4 Previous work
Confirmation that the features of Figure 2 are not numerical artefacts comes reports of from similar behaviour in Verrier & Evans (2009). In that work the authors report discovering a counterplay between the Kozai mechanism and a new circumbinary libration whilst modelling orbits within the double binary system HD98800.
Following on from Verrier & Evans (2009) is the excellent analytic paper of Farago & Laskar (2010). In this article the authors consider the elliptically restricted threebody problem and take advantage of the assumption that, in the circumbinary case, the displacement of the third body is greater than the relative separation of the binary .
Under this approximation Farago & Laskar expand the three body Hamiltonian to second order in and then average over the orbit of the binary and the third body to obtain a timeaveraged quadrupolar Hamiltonian.
This Hamiltonian promises to be very accurate in the regime , as higher order terms in will tend to zero at a faster rate than those of second order. This model should not be so good for orbits closer to the binary system, where higher order terms will play a more substantial role. In the following sections we investigate the accuracy of Farago & Laskar’s model by testing the predictions that it makes, and the limits at which their quadrupolar approximation becomes insufficient.
3.5 Separatrix and critical angle
The separatrix (Figure 3: black) is the boundary in phasespace between different modes of behaviour. In our circumbinary surface of section the separatrix takes the form of a triple figureofeight, or two circles intersecting at (, ) and (, ), dividing the regions of behaviour outlined in § 3.2 above. The separatrix which we plot in Figure 3 is actually the path of a rare test particle which, due to our stepping integrator, sampled more than one region of behaviour.
In Section 3.3 we mentioned that the geometry of the surface of section is predominantly dependent on binary eccentricity, and here we quantify this. We note that each point on the surface of section (Figure 2) defines a unique path, and that every path intersects the vertical axis (). We define a critical angle as the inclination at which the separatrix crosses the positive vertical axis () in the region .
Subsequent to the discussion of symmetry in § 3.3.2 it follows that the three other intersections of the separatrix with the vertical axis are reflections of the critical angle defined in the region (, ).
Measuring the critical angle
We have run an additional suite of simulations to find and extract the critical angle as a function of binary eccentricity. These simulations were run with one shell of test particles at radius , longitude of the ascending node and with a high resolution in inclination, intervals of , at various values of binary eccentricity.
For each orbit sampled we visually identify the behaviour (§ 3) as a function of inclination and binary eccentricity to find upper and lower boundaries on the separatrix. The results are plotted in Figure 4 preserving the colour scheme of Figure 2.
The green, red and blue points in Figure 4 represent the behaviour of test particles at these locations. We only plot the sampled points that lie either side of the separatrix. We extract the simplest polynomial fit to accurately describe the critical angle () as a function of binary eccentricity with the constraint that the fit must pass between every pair of greenred points. This is given by:
(1) 
where
In all cases the angle of the separatrix in the region () is at inclination , as it should be by symmetry arguments.
Comparison with Farago & Laskar (2010)
The timeaveraged quadrupolar model of Farago & Laskar (2010) predicts the critical angle (their equation 2.34) and hence the location of the separatrix as
(2) 
Since we measure the critical angle at a radius of from the binary we expect the quadrupolar approximation to hold, and indeed we find an excellent agreement between our measurements and Equation 2. The predicted location of the separatrix lies between every pair of our experimental limiting points.
We plot our data, fit and the Farago & Laskar prediction together in Figure 5. There is an almost exact agreement between model and data at this radius. The largest divergence is at , which is an unphysical limit as the binary itself becomes unbound.
3.6 Constant of motion
Verrier & Evans (2009) proposed an integral of motion (see their equation 3) for the libration islands to be the component of the specific angular momentum of an orbit along the line of apses of the internal binary, which may be expressed as
(3) 
This makes intuitive sense as the libration islands are centred at (, ), i.e. the points at which a test particle’s angular momentum is exactly parallel or antiparallel to the line of apses of the internal binary (see Figure 3). So for small deviations from these central points we find that is conserved.
Unfortunately this model breaks down as we move out from the centre of libration. In Figure 6 (upper panel) we show the Farago & Laskar constant of motion (Eq 3) over the course of our integration for example test particles in proximity to the centre of the island of libration (, ) from our simulation of binary eccentricity 0.6 and mass fraction 0.5. We observe that this suggested constant of motion becomes insufficient for test particles of lower inclination, which sample phase space away from the centre of the libration.
A more promising integral of motion is given by the timeaveraged quadrupolar model of Farago & Laskar (2010), their equation 2.20, which we translate to orbital elements as
(4) 
We see that this equation reduces to the square of Equation 3 for values of (, ) close to the centre of libration. But this model correctly predicts the paths of orbits across the entirety of the surface of section and for values of binary eccentricity . We plot this constant of motion in the lower panel of Figure 6 for the same test particles as above.
We use the constant of motion predicted by the timeaveraged quadrupolar model of Farago & Laskar (2010) as a test of the limits of the quadrupolar approximation. For each simulated test particle we have snapshots over the course of integration. We therefore calculate an instantaneous for each test particle at each snapshot and ask the question: ‘how constant is the constant of motion?’
For example — in our shell of test particles at radius we find that a typical test particle experiences variation in of . This is measured by taking the standard deviation of the instantaneous measurements of over the course of a simulation. Since is of order unity, this equates to a typical error of 0.0014% at these large radii.
We now investigate how well the quadrupolar approximation holds as we consider orbits closer to the binary. To do so we create histograms of for each shell of test particles that we sample (see Table 2). In Figure 7 we show histograms of for radii 3, 4, 5 and 6. We also colour the histograms according to the simulated binary eccentricity .
As we consider orbits closer to the binary (), we see that the quadrupolar approximation begins to break down. At 6 most orbits are accurate to the quadrupolar model to better than 0.5%, but as we move in to 3 some orbits experience deviation of . We observe similar small perturbations in a test particle’s semimajor axis and its initially zero eccentricity .
We also note that orbits around binary systems of higher eccentricity (Figure 7: purple) deviate significantly more than those around circular binaries (Figure 7: light blue), which indicates that the higher order terms of the Farago & Laskar Hamiltonian have a greater dependence on binary eccentricity.
Significance of
In Figure 8 we plot four typical examples of the variation in over the course of integration. These are selected from a random sample of 1000 such plots of test particles from the 3 radius bin (Figure 7, upper panel).
Numerical noise from a stepping integrator would in principle accumulate over time. In Figure 8 we show that the variation in is present from the start and constant in magnitude over the duration of a simulation. Values of are therefore not biased by an accumulation of numerical noise towards the end of the simulation.
3.7 Period of precession
From our finely timesampled data (Table 2) we extract a period for the precession of each stable test particle. These range from as little as (or ) to greater than the simulation duration, but good fits are obtained for the vast majority.
In Figure 9 we plot test particle traces on the surface of section for our simulation of binary eccentricity and binary mass fraction , as a function of test particle radius , colouring the traces by precession period.
We observe that the period of precession is correlated with a test particle’s orbital radius and proximity to a separatrix. Orbits which appear to be missing from Figure 9 are unstable, and hence are not plotted. We explore the topic of stability in much greater detail in the second part of the paper (§ 4).
As discussed above and demonstrated by Figure 9 the closer a circumbinary test particle orbits to the binary system the shorter is its period of precession . We find that a powerlaw provides a very good fit to as a function of distance from the binary, and so we fit to all orbits sampled, where is a free parameter. An example fit is shown in Figure 10. The data point from our simulation at radius provides a very good constraint.
Our data show a tight clustering about , as plotted in Figure 11. This is not within the range given by Verrier & Evans (2009) of , but does agree with the timeaveraged quadrupolar model of Farago & Laskar which predicts of exactly 3.5.
But on closer inspection of Figure 11 we notice that appears to consist of two populations. One of these populations is comprised of the prograde and retrograde orbits, which cluster around , and the other corresponds to the polar orbits, which cluster around .
An analytic expression for period of precession
The timeaveraged quadrupolar model of Farago & Laskar (2010) makes the following prediction (their equation 2.32):
(5) 
where as defined in Eq 4,
and may be defined in terms of the complete elliptical integral of the first kind, , as
where
We calculate the expected period (Eq 5) for each test particle that we sample, and in Figure 12 we plot a straight onetoone comparison between the predicted and measured periods of precession.
Figure 12 shows two populations — the upper right group consists of test particles from our radius simulations, hence the larger periods, whereas the lower left group represents every test particle from the main body of our simulations (1.5, summarised in Table 2).
The scatter in period at the upper end of the groups of Figure 12 is due to these periods being significantly greater than the duration of the simulation, and hence badly fitted. The scatter at lower periods () is due to deviations of the experimental values from the quadrupolar modal at low radii. But within the well behaved boxed region of Figure 12:
The results of our numerical investigation into the dynamics of circumbinary orbits in the elliptically restricted three body problem have correlated well with the analytic work of Farago & Laskar (2010). We have also explored the limits of their model due to the quadrupolar approximation.
4 Stability of circumbinary orbits
Figure 9 shows test particle traces on the surface of section for our simulation of binary eccentricity and binary mass fraction , as a function of test particle orbital radius from the centre of mass of the binary. Each radius bin is equally sampled but we only plot orbits which remained stable throughout the simulation. It can be seen that at a radius of only very few of the initial test particles are actually stable, but at greater distances from the binary the surface of section fills in.
Figure 9 shows that in the plotted simulation (, ) the closest stable orbits to the binary are those at the centres of the islands of libration. These orbits are stable at a radius of , with regions closer to the separatrix becoming stable out to . As we move further away from the binary the retrograde orbits begin to acquire stability for radii , before the prograde orbits finally become stable at .
The timeaveraged quadrupolar model of Farago & Laskar (2010) predicts the dynamics of circumbinary orbits but can make no attempt to discern whether these orbits are viable. In the following sections we reveal characteristics which are due to resonances between the binary and test particle orbital periods which the Farago & Laskar model cannot explain due to the timeaveraging carried out in its derivation.
4.1 Suite of simulations
We ran a suite of simulations to investigate the stability of circumbinary orbits. As discussed in § 3.5 the dynamics of the circumbinary phase space are such that every orbit crosses the axis and there exists a symmetry which reflects onto . We may therefore narrow down the region of phase space which we sample to only one value of the longitude of the ascending node, .
We run these simulations to 50,000 binary orbital periods, sampling phase space to a good density in orbital radius and inclination , as laid out in Table 3. We sample across binary eccentricity and mass fraction as in the above sections (see Table 1).
Orbital Element  min  max  

semimajor axis  0.05  
inclination  
longt. of the asc. node  
true anomaly  
sim. length and snapshot 
= binary semimajor axis and = binary orbital period.
4.2 A measure of stability
As discussed in § 2.2, each test particle is monitored for instability. Unstable orbits are identified and removed during integration where a test particle is perturbed sufficiently from its initial orbit to approach either star, or if it evolves onto an unbound trajectory (). Postsimulation stability criteria are applied to identify and reject test particles which do not quite reach escape velocity.
Whereas in section 3 we were concerned with the orbits which survived the simulation, here we are more interested in those which don’t. In Figure 13 we plot a histogram of the escape times at which particles are rejected from the simulation. The vast majority of unstable particles are caught at the start of the simulation — of the approximately half a million unstable test particles, over half of these are caught within the first 1000 binary orbital periods, with the distribution tailing off steeply even in logspace.
Since we sample each orbit from multiple initial values of the true anomaly (Table 1) we measure an orbit’s longterm stability by the fraction of initial test particles which survive the simulated duration of .
In Figure 14 we show a density plot of this measure of stability across our entire parameter space. The major axes of this figure correspond to the simulation parameters of binary eccentricity and mass fraction, whilst the minor axes correspond to orbital radius and inclination. This is a density plot where each pixel corresponds to a sampled orbit, and the transparency to our measure of stability — a darker colour indicates a more stable orbit. We preserve the colour scheme of Figure 2.
We draw the reader’s attention to the following features of Figure 14, as revealed by our exquisitely detailed simulations:

First, and very broadly speaking, orbits are more stable at lower binary eccentricity .

With equal generality, retrograde orbits (blue) appear to be the most stable, followed by librating orbits (red), and finally prograde orbits (green). The difference in radius of the innermost stable orbit across inclination can be as large as .

There are vertical striations of instability, most noticeable in the higher eccentricity simulations , and predominantly at inclinations of and . We hypothesize that these regions of instability are due to orbital resonances between a test particle and the binary.

We note very thin horizontal pinnacles of instability in nonlibrating orbits . These pinnacles are located at inclinations and , and extend up to into otherwise stable phase space.

More central to the pinnacles discussed above are wider peninsulas of instability in the librating region. These appear symmetrically either side of in the librating region for simulations of and converge upon each other as .

The horizontal pinnacles and peninsulas are a function of binary mass fraction . These features do not appear in the simulations, and are magnified towards increasingly extreme values of .
4.3 Previous work
There exist numerous papers in the literature in which authors investigate longterm orbital stability within the coplanar circular restricted three body problem, both numerically and analytically, and with emphasis on circumstellar and circumbinary orbits.
With improvements in computation power authors have relaxed the circular constraint in the problem to investigate eccentric binary systems (Dvorak et al., 1989; Holman & Wiegert, 1999; Musielak et al., 2005). But only recently have we had the computation power to relax the coplanar constraint on the problem to investigate inclined orbits.
PilatLohinger
et al. (2003) performed three dimensional numerical experiments to determine inclined stability but they did not explore the complex libration structure of the phase space
This paper presents the first dynamicaware analysis of the stability of inclined circumbinary orbits throughout binary mass fraction — eccentricity parameter space.
4.4 Escape time
In a companion figure to the stability plot of Figure 14 we show the escape time of each unstable orbit in Figure 15. Here we use the same axes as Figure 14 — with major axes corresponding to the simulation parameters of binary eccentricity and mass fraction, and minor axes corresponding to orbital radius and inclination. This is a density plot where each pixel corresponds to an orbit sampled, and the transparency to the inverse of escape time — a darker colour indicates a longer surviving orbit. We again preserve the colour scheme of Figure 2. Since we consider multiple test particles per sampled orbit we take an average for the escape time.
In Figure 15 we find matching features to those in Figure 14, as described above in § 4.2. But here we can also observe how longlived the unstable orbits are. For example, in the low binary eccentricity simulations we find that orbits are either very quickly unstable, or definitely stable. But for the simulations of higher binary eccentricity and more imbalanced binary mass fraction , the most striking features are longlived. The test particles which sample phase space at these extreme points are almost stable, and survive for times of order the simulated duration (Table 3).
Test particles on unstable orbits further inside the unstable region (closer to the binary) are quickly accreted onto the binary stars and/or ejected from the system. It is within this region of instability that inflows and outflows are likely to be important.
4.5 Further considerations
The phase space underlying our study of stability within the threebody problem is arguably significantly chaotic. We believe that the features which we report are real and significant to the consideration of matter around binaries, but we advise caution before cranking up the phasespace resolution. We should consider the effect of perturbations away from this idealised model, such as the volume of the bodies considered and the feedback effect on the internal binary of a third body of nonnegligible mass.
5 Conclusions
Our simulations show that inclined circumbinary orbits in the ellipticallyrestricted threebody problem demonstrate three distinct families of behaviour: closetocoplanar prograde () and retrograde () orbits precess in the longitude of the ascending node, whilst closetopolar orbits ( and ) have their longitude of the ascending node and inclination coupled to precess about the centre of an island of libration.
We have extracted the critical angle of the separatrix along the critical axis between these regions of behaviour as a function of binary eccentricity.
We have shown that the analytic timeaveraged quadrupolar model of Farago & Laskar (2010) provides an excellent description of the behaviours of circumbinary orbits at radii . We have also shown that their model becomes inaccurate to greater than 1% at orbital radii , and especially in cases of high binary eccentricity.
With the first 3D dynamicaware analysis of the stability of circumbinary orbits we have discovered that these orbits are surprisingly stable throughout binary mass fraction — eccentricity parameter space. Our detailed simulations have put numerical limits on this stability and revealed complex structure in orbital radius — inclination space.
Our work shows that circumbinary phasespace is rich and dynamic, full of remarkable and stable orbits which do not behave simply. We should not presume any given binary system to lack a circumbinary component unless otherwise demonstrated. Such a component may be a source of obscuration, emission, inflow or outflow.
Acknowledgments
SD thanks STFC for a studentship and KMB thanks the Royal Society for a University Research Fellowship. We thank John Magorrian for useful discussions and our referee for his careful reading and useful comments towards this manuscript.
Footnotes
 i.e. the binary system’s semimajor axis
 ‘Ptype’ (Dvorak et al., 1989)
 An orbit about one star of a binary system. ‘Stype’ (Dvorak et al., 1989).
 yet remaining in the binary orbital plane
 The inclined simulations of PilatLohinger et al. contain only test particles of initial longitude of the ascending node (PilatLohinger et al., private communication). As such they did not sample any librating orbits, which do not intersect the or axes.
References
 Artymowicz P., Lubow S. H., 1996, \apjl, 467, L77+
 Beuermann K., Buhlmann J., Diese J., Dreizler S., Hessman F. V., Husser T.O., Miller G. F., Nickol N., Pons R., Ruhr D., Schmülling H., Schwope A. D., Sorge T., Ulrichs L., Winget D. E., Winget K. I., 2011, \aap, 526, A53+
 Blundell K. M., Bowler M. G., Schmidtobreick L., 2008, \apjl, 678, L47
 Doolin S., Blundell K. M., 2009, \apjl, 698, L23
 Dvorak R., Froeschle C., Froeschle C., 1989, A&A, 226, 335
 Farago F., Laskar J., 2010, \mnras, 401, 1189
 Green R. M., 1985, Spherical Astronomy. Cambridge University Press
 Holman M. J., Wiegert P. A., 1999, \apjl, 117, 621
 Kozai Y., 1962, \aj, 67, 591
 Lee J. W., Kim S., Kim C., Koch R. H., Lee C., Kim H., Park J., 2009, \apj, 137, 3181
 Musielak Z. E., Cuntz M., Marshall E. A., Stuit T. D., 2005, A&A 434, 355364
 Perez M. S., Blundell K. M., 2010, \mnras, 408, 2
 PilatLohinger E., Funk B., Dvorak R., 2003, A&A, 400, 1085
 Qian S.B., Liu L., Liao W.P., Li L.J., Zhu L.Y., Dai Z.B., He J.J., Zhao E.G., Zhang J., Li K., 2011, \mnras, pp L241+
 Regös E., Bailey V. C., Mardling R., 2005, \mnras, 358, 544
 Schwarz R., Haghighipour N., Eggl S., PilatLohinger E., Funk B., 2011, ArXiv eprints
 Verrier P. E., Evans N. W., 2009, \mnras, 394, 1721