Orbital dynamics of circumbinary planets

ABSTRACT

We investigate the dynamics of a non-zero mass, circular orbit planet around an eccentric orbit binary for various values of the binary eccentricity, binary mass fraction, planet mass, and planet semimajor axis by means of numerical simulations. Previous studies investigated the secular dynamics mainly by approximate analytic methods. In the stationary inclination state, the planet and binary precess together with no change in relative tilt. For both prograde and retrograde planetary orbits, we explore the conditions for planetary orbital libration versus circulation and the conditions for stationary inclination. As was predicted by analytic models, for sufficiently high initial inclination, a prograde planet’s orbit librates about the stationary tilted state. For a fixed binary eccentricity, the stationary angle is a monotonically decreasing function of the ratio of the planet-to-binary angular momentum j. The larger j, the stronger the evolutionary changes in the binary eccentricity and inclination. We also calculate the critical tilt angle that separates the circulating from the librating orbits for both prograde and retrograde planet orbits. The properties of the librating orbits and stationary angles are quite different for prograde versus retrograde orbits. The results of the numerical simulations are in very good quantitative agreement with the analytic models. Our results have implications for circumbinary planet formation and evolution.

methods: analytical, celestial mechanics, binaries: general

1 INTRODUCTION

Binary stars form in turbulent molecular clouds (McKee & Ostriker 2007) where the accretion process during star formation may be chaotic (Bate, Bonnell & Bromm 2003; Bate 2018). Therefore circumbinary discs likely form misaligned with respect to the binary orbital plane. Observationally, misaligned circumbinary discs appear common for wider binaries (e.g. Chiang & Murray-Clay 2004; Winn et al. 2004; Kennedy et al. 2012, 2019; Brinch et al. 2016). Observations suggest that there is a break between alignment and misalignment of circumbinary discs at orbital periods of about |$30 \, \rm d$| (Czekala et al. 2019). For short-period binaries, period |$\lt 20\, \rm d$|⁠, 68 per cent have aligned discs (within 3°). Giant planets form while the gas disc is still present (e.g. Lagrange et al. 2010) and thus it is likely that planets may form on misaligned orbits. A giant planet in a misaligned disc in a binary system may not remain coplanar to the disc (Picogna & Marzari 2015; Lubow & Martin 2016; Martin et al. 2016; Franchini, Martin & Lubow 2019a).

In this paper, we concentrate on the properties of circumbinary planets involving eccentric orbit binaries. The Kepler mission has so far detected 10 circumbinary planets. Among these, at least two of them have been detected around eccentric binaries. Kepler-34b with a mass of |$0.22\, M_{\rm J}$| orbits an eclipsing binary star system (Kepler-34) that has an orbital eccentricity of 0.52 (Welsh et al. 2012; Kley & Haghighipour 2015). The moderately eccentric (e = 0.26) binary KIC 5095269 hosts a circumbinary planet KIC 5095269b (Getley et al. 2017).

All the circumbinary planets detected thus far are nearly coplanar with the binary orbital plane (e.g. Li, Holman & Tao 2016). However, this is likely a selection effect due to the small orbital period of the Kepler binaries. Longer orbital period binaries are then expected to host planets with a wide range of inclinations. Planets with large misalignments are much more difficult to detect than those that orbit in the binary orbital plane because their transits are much more rare. However, other detection methods may be possible. For example, polar planets (those in an orbit close to perpendicular to the binary orbital plane) may be distinguished from coplanar planets through eclipse timing variations (Zhang & Fabrycky 2019).

For a misaligned test (massless) particle circumbinary orbit about a circular orbit binary, its nodal precession occurs around the binary angular momentum vector, but may be either prograde or retrograde depending upon the initial particle inclination. The nodal precession about a circular orbit binary is always circulating. That is, the longitude of the ascending node fully circulates over 360°, since the nodal precession rate does not change sign. However for a binary with non-zero eccentricity, a circumbinary test particle orbit with a sufficiently large inclination may undergo libration. That is, the longitude of the ascending node covers a limited range of angles less than 360° and the nodal precession rate changes sign. In the test particle case, the angular momentum vector of the test particle librates about the binary eccentricity vector (or binary semimajor axis) and undergoes tilt oscillations (Verrier & Evans 2009; Farago & Laskar 2010; Doolin & Blundell 2011; Naoz et al. 2017; de Elía et al. 2019). The minimum inclination required for libration decreases with increasing binary eccentricity. This means that a test particle orbit with even a small inclination can librate around a highly eccentric binary.

Misaligned low-mass/angular momentum discs in and around binaries can undergo similar precession as a test particle (Larwood et al. 1996). Consequently, following the behaviour of test particles, a sufficiently misaligned low-mass disc around an eccentric binary can precess around the eccentricity vector, rather than the binary angular momentum and undergo tilt oscillations (Aly et al. 2015; Martin & Lubow 2017, 2018; Lubow & Martin 2018; Zanazzi & Lai 2018; Franchini, Lubow & Martin 2019b). Dissipation within the disc leads to the eventual alignment, either coplanar or polar aligned with respect to the binary orbital plane. In the polar aligned state, the low angular momentum circumbinary disc lies perpendicular to the binary orbital plane and its angular momentum vector is along the binary eccentricity vector. The polar disc as well as the binary do not undergo nodal precession. In the eccentric orbit binary, a disc that is evolving towards coplanar alignment undergoes tilt oscillations as it does so (Smallwood et al. 2019).

The mass of the disc has a significant effect on the polar alignment. A disc with mass is expected to evolve to a generalized polar state in which the inclination of the disc relative to the binary is stationary in a frame that precesses with the binary. A simplified model for the disc with mass involves the orbit of a particle with mass. The particle orbit then corresponds to a ring or narrow disc. In this stationary state, the disc inclination is less than 90° for a disc in a prograde orbit (Zanazzi & Lai 2018; Martin & Lubow 2019).

In this work, we extend the three-body numerical simulations of a circular orbit, test (massless) circumbinary particle performed by Doolin & Blundell (2011) by allowing the particle (planet) to have non-zero mass. There has been some exploration into the stability of low-mass polar particles (Cuello & Giuppone 2019). In this case, the binary feels the gravitational force of the planet, which causes the binary orbit to be modified. The eccentricity vector of the binary precesses as a result of this interaction. The binary orbit tilt and the magnitude of its eccentricity oscillate in this case. We consider the third body (planet) to lie on a circular circumbinary orbit. The dynamics of a circular orbit planet are similar to those of a narrow ring that in turn may be indicative of a more extended disc with mass. Therefore, the results can also have implications for the evolution of a circumbinary disc with non-zero mass. Our simulations are scale free and can therefore be applied to all scales.

In Section 2, we describe the initial conditions for our three-body simulations and explore the properties of circumbinary orbits with a planet, represented as a particle with non-zero mass. In Section 3, we determine the critical angles between the librating and circulating orbits and compare them to the analytic solutions found by Martin & Lubow (2019) that is based on earlier work by Farago & Laskar (2010) who used the quadrupole approximation. These orbits have been studied in the octupole approximation only for test particles (Naoz et al. 2017; Vinson & Chiang 2018). We also determine stationary inclinations, the inclination at which the binary and the third body orbit precess at the same rate with no change in relative tilt. In Section 4, we present our discussion and conclusions.

2 THREE-BODY SIMULATIONS

In this section, we describe and present results of our three-body simulations. We explore the effects of the binary eccentricity and the binary mass ratio on the different orbit types, circulating and librating, for planets with varying angular momentum.

2.1 Three-body simulation set-up

To study the evolution of a third body orbiting around an eccentric binary star system, we use the N-body simulation package, rebound. We use a WHfast integrator which is a second-order symplectic Wisdom Holman integrator with 11th-order symplectic correctors (Rein & Tamayo 2015). We solve the gravitational equations for the three bodies in the frame of the centre of mass of the three-body system. The central binary has components of mass m₁ and m₂ with total mass m_b = m₁ + m₂ and the mass fraction of the binary is f_b = m₂/m_b. The orbit has semimajor axis a_b, the magnitude of the eccentricity of the binary is e_b, and the orbital period of the binary is T_b. The Keplerian orbit of the planet with mass m_p around the centre of mass of the binary is defined initially by six orbital elements, its semimajor axis a, inclination i, eccentricity e, longitude of the ascending node ϕ, argument of periapsis ω, and true anomaly ν. Since the planet orbit is initially circular, we set as initial conditions e = 0 and ω = 0. We take ν = 0 and ϕ = 90° initially in our suites of simulations. We vary the planet mass, initial inclination, and semimajor axis of its orbit. Note that the binary orbit is not fixed since the binary feels the gravity of the massive third body.

In order to plot the results, we work in a frame defined by the instantaneous values of the eccentricity and angular momentum vectors of the binary, |$\boldsymbol {e}_{\rm b}$| and |$\boldsymbol {l}_{\rm b}$|⁠, respectively. The frame defined by the binary has the three axes |$\boldsymbol {e}_{\rm b}$|⁠, |$\boldsymbol {l}_{\rm b}\times \boldsymbol {e}_{\rm b}$|⁠, and |$\boldsymbol {l}_{\rm b}$|⁠. Denoting the planet angular momentum as |$\boldsymbol {l}_{\rm p}$|⁠, we determine the inclination of its orbital plane relative to the binary through

$$\begin{eqnarray*} i =\cos ^{-1}(\boldsymbol {\hat{l}}_{\rm b}\cdot \boldsymbol {\hat{l}}_{\rm p}), \end{eqnarray*}$$

(1)

where |$\boldsymbol {\hat{l}}_{\rm b}$| is a unit vector in the direction of the angular momentum of the binary and |$\boldsymbol {\hat{l}}_{\rm p}$| is a unit vector in the direction of the angular momentum of the particle. The inclination of the binary relative to the total angular momentum |$\boldsymbol {l}$| is

$$\begin{eqnarray*} i_{\rm b} =\cos ^{-1}(\boldsymbol {\hat{l}}\cdot \boldsymbol {\hat{l}}_{\rm b}), \end{eqnarray*}$$

(2)

where |$\boldsymbol {\hat{l}}$| is a unit vector in the direction of the total angular momentum (⁠|$\boldsymbol {l}=\boldsymbol {l}_{\rm b}+\boldsymbol {l}_{\rm p}$|⁠). Similarly, we determine the phase angle of the particle in the same frame of reference through

$$\begin{eqnarray*} \phi =\tan ^{-1}\left(\frac{\boldsymbol {l}_{\rm p}\cdot (\boldsymbol {l}_{\rm b}\times \boldsymbol {e}_{\rm b})}{\boldsymbol {l}_{\rm p}\cdot \boldsymbol {e}_{\rm b}} \right). \end{eqnarray*}$$

(3)

In the next three subsections, we vary masses and orbital properties of the binary stars and planets and describe the binary and planet orbital evolution. We first ran some test particle simulations with the planet mass set to zero in order to verify that our results are in agreement with the results presented in fig. 2 in Doolin & Blundell (2011). We then performed a set of simulations for various values of the mass of the planet from |$m_{\rm p}=0.001\, m_{\rm b}$| to |$m_{\rm p}=0.01\, m_{\rm b}$| and its orbital radius from |$r=5\, a_{\rm b}$| to |$r=20\, a_{\rm b}$|⁠. We also explored the effect of different initial binary eccentricities and mass ratios. We also consider some simulations with higher angular momentum of the third body. In Table 1, we list the parameters for each model.

Table 1.

Parameters of the simulations. The first column contains the name of the model, the second and third columns indicate the binary mass fraction and initial eccentricity. The fourth and fifth columns represent the mass of the planet in units of m_b and the distance of the planet with respect to the centre of the mass in units of a_b, respectively.

Model	fb	eb	mp/mb	r/ab
A1	0.5	0.2	0.001	5
A2	0.5	0.5	0.001	5
A3	0.5	0.8	0.001	5
B1	0.1	0.2	0.001	5
B2	0.1	0.5	0.001	5
B3	0.1	0.8	0.001	5
C1	0.5	0.2	0.01	5
C2	0.5	0.5	0.01	5
C3	0.5	0.8	0.01	5
D1	0.1	0.2	0.01	5
D2	0.1	0.5	0.01	5
D3	0.1	0.8	0.01	5
E1	0.5	0.2	0.01	20
E2	0.5	0.5	0.01	20
E3	0.5	0.8	0.01	20
F1	0.1	0.2	0.01	20
F2	0.1	0.5	0.01	20
F3	0.1	0.8	0.01	20
G1	0.5	0.2	0.056	20
G2	0.5	0.5	0.050	20
G3	0.5	0.8	0.034	20
H1	0.5	0.2	0.116	20
H2	0.5	0.5	0.102	20
H3	0.5	0.8	0.070	20

2.2 Low-mass planet at small orbital radius

Fig. 1 shows the results of our three-body simulations for a low-mass planet with |$m_{\rm p}=0.001\, m_{\rm b}$| on an orbit with semimajor axis |$r=5\, a_{\rm b}$| for varying e_b and f_b. Each line in each plot corresponds to a planet orbit with different initial inclination to the binary. The first and third columns show the icos ϕ–isin ϕ plane while the second and fourth columns show the corresponding e_bcos ϕ–e_bsin ϕ plane. The black stars represent the initial position of the planet.

Figure 1.

The icos ϕ–isin ϕ plane (first and third columns) and e_bcos ϕ–e_bsin ϕ plane (second and fourth columns) for orbits with different values of initial inclination and longitude of the ascending node. The planet has mass |$m_{\rm p} = 0.001\, m_{\rm b}$|⁠, and orbital radius |$5\, a_{\rm b}$|⁠. The binary eccentricity is |$e_{\rm b}=0.2,\, 0.5\, ,0.8$| in the upper, middle, and lower panels, respectively. The mass fraction of the binary is f_b = 0.5 in the first and second column and f_b = 0.1 in the third and fourth column. The green and blue lines represent prograde and retrograde circulating orbits, respectively. The red lines represent librating orbits with initial inclination i < i_s while the cyan lines represent librating orbits with initial inclination i > i_s. The black stars mark the initial positions of the planet with i ranging from 10° to 180°. We removed unstable orbits.

The green lines represent prograde (relative to the binary) circulating orbits where the planet displays clockwise precession in the longitude of the ascending node ϕ. The blue lines correspond to retrograde circulating orbits where the planet displays counterclockwise precession in ϕ. The red and cyan lines identify librating orbits. The inclination at the centre of these orbits is the stationary inclination, i_s. In this low-mass planet case, the centres are at i = i_s ≈ 90° and ϕ = ±90°. The red lines have initial inclination i < i_s while the cyan lines have initially i > i_s. These librating orbits display counterclockwise precession in ϕ.

The icos ϕ–isin ϕ phase plots are very similar to those in the test particle case considered by Doolin & Blundell (2011). The eccentricity of the binary does not vary during the precession of a test particle because the particle does not have angular momentum to exchange with the binary system. The e_bcos ϕ–e_bsin ϕ panels for these models show the curves are sometimes circulating while slightly non-circular and sometimes librating due to the eccentricity oscillations associated with the relatively small but non-zero particle (planet) mass. If instead the object orbiting the binary is a planet (a non-zero mass object), e_b initially increases or decreases depending on whether the initial angle of the planet’s orbit is greater or smaller than the stationary inclination. As we shall see, the larger the planet angular momentum, the larger the oscillations of the binary orbit.

In order to investigate the effect of the initial binary eccentricity, we ran simulations with initial e_b  = 0.2, 0.5, and 0.8. Comparing the three rows in Fig. 1, we see that higher initial eccentricities correspond to larger libration islands.

The comparison between Models A and B in Fig. 1 shows the effect of decreasing the binary mass fraction from f_b = 0.5 to f_b = 0.1. The size and shape of the orbits in the icos ϕ–isin ϕ plane do not change significantly with mass fraction, but the variations of e_b in the libration region are larger for smaller binary mass fraction. We find that for lower binary mass fraction, the orbits near the libration centre (with inclination i = i_s and phase angle ϕ = ±90°) become divergent for moderate to high binary eccentricities (e_b  = 0.5, 0.8) and therefore might be unstable. This issue is beyond the scope of this work. We will explore the stability of orbits close to misaligned binaries in a future publication.

In the top row of Fig. 2, we consider in more detail the evolution of the orbits in time for models with initial binary eccentricity of 0.5 and binary mass fraction f_b = 0.5 (top left, model A2) and f_b = 0.1 (top right, model B2). We show e_b, i_b, and i as a function of time for different values of the initial planet inclination. The colour of the lines corresponds to the orbit type in Fig. 1. Comparing the left and right plots, we see that there is less evolution of the binary in the equal mass binary compared to the lower binary mass fraction. This is because with lower binary mass fraction, the planet has more angular momentum compared to the binary and thus has a stronger effect on the binary. The time-scale for the oscillations is longer for smaller binary mass fraction.

Figure 2.

Time evolution of the binary eccentricity, e_b, the inclination of the binary i_b with respect to the vector of the total angular momentum, and the inclination of the planet with respect to the vector of the binary angular momentum,i, for Models A2, B2, C2 D2, E2, and F2. Each panel contains one line for each different type of orbit with the difference being the initial inclination, which is shown in the bottom panel. The colours correspond to the orbit type described in Fig. 1.

2.3 High-mass planet at small orbital radius

Fig. 3 shows the results in the icos ϕ–isin ϕ and e_bcos ϕ–e_bsin ϕ plane obtained with a higher mass planet |$m_{\rm p} = 0.01\, m_{\rm b}$| orbiting an eccentric binary at |$r=5\, a_{\rm b}$|⁠. The line colours correspond to the same type of orbits as in Fig. 1. Comparing Fig. 3 with Fig. 1, we see that the region of prograde circulating orbits is larger for a higher mass planet. The differences between the red and cyan lines becomes more prominent. The inclination at the centre of the libration, i_s, decreases. For example, for model C2 with e_b = 0.5 and f_b = 0.5 it is 80°. Because of this decrease, there is a narrower range of inclinations for which red orbits exist, those with initial inclination i < i_s, and a wider range of inclinations for which cyan orbits exist with i > i_s initially. We discuss the value of the stationary inclination further in Section 3.

Figure 3.

Same as Fig. 1 except |$m_{\rm p} = 0.01\, m_{\rm b}$|⁠.

The higher planet mass causes the binary eccentricity to vary significantly not only in the librating solutions, but also in the circulating orbits of the planet. The binary eccentricity in the libration islands for Model C2 in Fig. 3 starts at e_b = 0.5 and reaches values as large as e_b ≃ 0.7.

By comparing Models C and D in Fig. 3, we can see the effect of changing the binary mass fraction. The cyan libration islands are significantly larger and there are fewer red orbits in the simulation with f_b  = 0.1 because i_s is smaller than in the equal mass binary case. Therefore, decreasing the binary mass fraction results in a lower stationary inclination. We find particle orbital instability close to i = i_s and ϕ = 90° for Model D1. There are no stable orbits in the region close to the stationary inclination.

The variations in e_b are larger for the smaller binary mass fraction, as seen in the e_bcos ϕ–e_bsin ϕ plot. In particular, e_b in Model D3 becomes very close to 1 during libration.

The middle panels of Fig. 2 shows the binary eccentricity, inclination, and planet inclination evolution with time for Models C2 and D2. We see again that decreasing the binary mass fraction leads to larger amplitude oscillations in the binary eccentricity for librating orbits starting above the critical angle (cyan line). Model D2 has larger variations of i_b and i during the libration of the planet because the system has a larger planet-to-binary angular momentum ratio and the secondary star (m₂ = 0.1m_b) interacts more strongly with the planet.

2.4 High-mass planet at large orbital radius

We now increase the planet semimajor axis to |$r=20\, a_{\rm b}$| and keep its mass at |$m_{\rm p}=0.01\, m_{\rm b}$|⁠. Comparing Model E in Fig. 4 with Model C in Fig. 3, we see that the libration islands become even larger when the planet orbits the binary with a larger semimajor axis. The eccentricity of the binary e_b can be excited to larger values during the planet’s libration because in this configuration the planet has more angular momentum to exchange with the binary system. The increase in the initial binary eccentricity e_b has the same effect as in previously described simulations, i.e. it increases the range of stable librating orbits.

Figure 4.

Same as Fig. 3 except the orbital radius is 20 a_b. The magenta lines show crescent orbits.

We find that there are no retrograde circulating solutions for the low-mass fraction case f_b = 0.1 (see the third column of Fig. 4, there are no blue lines).

The eccentricity of the binary can be excited to very large values (close to 1) during the planet’s precession in Models F1, F2, and F3. However, there is a new type of orbit that is different from those described in previous sections. These librating orbits, represented by magenta lines, only appear in the simulations with smaller binary mass fraction f_b = 0.1 that start with higher initial eccentricities |$e_{\rm b}=0.5,\, 0.8$|⁠. The magenta orbits have higher initial inclination than the librating cyan lines and display counterclockwise precession. In Model F2, these crescent shape orbits appear in the icos ϕ–isin ϕ plane for i > 150° while in Model F3, they appear also for lower inclinations i > 140°. These librating orbits are not nested within each other, as they are in the prograde case. These crescent orbits have been seen before in simulations, but not explored (e.g. Zanazzi & Lai 2018).

Appendix A of Martin & Lubow (2019) shows that non-coplanar stationary states for retrograde orbits with i_s < 180° do not exist below a critical value of the angular momentum ratio j_cr given by equation A3 in that paper

$$\begin{eqnarray*} j_{\rm cr}= \frac{1+4 e_{\rm b}^2}{2+3 e_{\rm b}^2}. \end{eqnarray*}$$

(4)

The lack of stationary states means that the crescent shaped librating orbits are always of non-zero extent in the icos ϕ–isin ϕ phase plane. This is understood by the lack of stationary retrograde inclinations for j < j_cr. For example, for Model F2 j = 0.58 < j_cr = 1.29 and for Model F3 j = 0.83 < j_cr = 0.91.

Appendix A of Martin & Lubow (2019) shows analytically that stationary coplanar retrograde i = 180° (and also prograde i = 0°) orbits should exist. However, we have not been able to find such an orbit numerically for Models F1, F2, and F3.

In the bottom panels of Fig. 2, we show the evolution of e_b, i_b, and i for Model E2 and F2. The time oscillations of e_b, i_b, and i have longer periods compared to the corresponding case of the same planet mass with smaller semimajor axis. We show the evolution of one of the crescent orbits, the magenta lines in Model F2. As e_b is very close to 1 during the precession, the vector of the binary angular momentum changes quickly, resulting in the narrow peaks in the i_b and i plots.

2.5 Systems with high angular momentum ratios

To investigate the retrograde librating orbits, we now consider simulations of an equal mass binary with higher angular momentum ratios j = 1 and j = 2, with binary eccentricities e_b  = 0.2, 0.5, and 0.8. The angular momentum ratios are larger than the critical required for a retrograde libration centre given in equation (4).

Fig. 5 shows the orbital evolution of simulations with j = 1 (left-hand panels) and j = 2 (right-hand panels). There are librating orbits in the phase diagrams in the first and third columns of Fig. 5 that surround the stationary points. with i_s < 90°. The fully retrograde librating orbits (i > 90° throughout the orbit) are seen as the crescent magenta orbits in the icos ϕ–isin ϕ phase plane. In the case of inclinations less than the critical inclination, the orbits decrease in extent in the phase plane with increasing initial inclination. They are at most only partially overlapping and are not fully nested. They reach zero extent at the stationary angle, i_s. For inclinations above the stationary angle, the orbits increase in extent in the phase plane. Above the critical angle, the binary eccentricity initially decreases while below the critical angle the binary eccentricity initially increases starting at ϕ = 90°, indicated by the stars in the phase planes. Unlike the prograde librating case, the retrograde librating orbits are not nested about a common centre that occurs at the stationary inclination.

Figure 5.

Same as Fig. 4 except that each model has different m_p to satisfy j = 1 (left) and j = 2 (right).

3 STATIONARY AND CRITICAL ANGLES

We compare the results of our numerical simulations to the analytic results presented in Martin & Lubow (2019) for the stationary inclination (inclination at the centre of the libration island), i_s, and the critical minimum inclination angle i_min that separates the prograde orbits from the librating orbits. We also calculate numerically the critical maximum inclination for librating orbits, i_max. The analytic results are based on secular equations with the quadrupole approximation for the binary potential (Farago & Laskar 2010). The equations are expected to break down for orbits that are close to the binary.

We note that the planet orbit remains nearly circular in our simulations, as is expected analytically for the simulations with f_b = 0.5 since the particle eccentricity is a constant of motion in the secular quadrupole approximation for the binary (Farago & Laskar 2010). We will discuss the particle eccentricities obtained from our simulations in this Section. Analytic calculations have been carried out to octupole order but only for circumbinary test particles (Naoz et al. 2017; Vinson & Chiang 2018).

3.1 Stationary inclination

The stationary inclination i_s depends only the eccentricity of the binary, e_b, and the ratio of the angular momentum of the particle to the angular momentum of the binary, j. The equation is given by

$$\begin{eqnarray*} \cos i_{\rm s} = \frac{-(1+4e_{\rm b}^2) \pm \sqrt{(1+4e_{\rm b}^2)^2+60(1-e_{\rm b})j^2}}{10j}\, \end{eqnarray*}$$

(5)

(see equation 17 in Martin & Lubow 2019), where we take the positive square root for prograde orbits (cos i_s > 0) and the negative square root for retrograde orbits (cos i_s < 0; see appendix A of Martin & Lubow 2019).

3.1.1 Prograde stationary inclination

The solid lines in Fig. 6 plot equation (5) in the prograde case for e_b = 0.2 (blue line), 0.5 (black line), and 0.8 (red line) as a function of the planet-to-binary angular momentum ratio j. For fixed e_b, the stationary inclination i_s decreases monotonically with increasing j and for fixed j it increases monotonically with increasing e_b. The stationary inclinations for the models in Table 1 were determined numerically from the simulations. The magenta dots in Fig. 6 correspond to the models with f_b = 0.5 and the cyan dots correspond to the models with f_b = 0.1. The results confirm the prediction that for fixed binary eccentricity the stationary inclination depends directly on j, independent of f_b. The simulations are in very good agreement with the analytically results. Thus, the higher order approximations such as octupole are not required to explain these results.

Figure 6.

Comparison of the analytic solution given by equation (5) in the prograde case (solid lines) with simulation results (dots) for the stationary tilt i_s of the planet relative to the binary as a function of planet-to-binary angular momentum ratio j. The solid curves have binary eccentricities e_b = 0.2 (blue line), 0.5 (black line), and 0.8 (red line). Cyan dots correspond to models with f_b = 0.1 and magenta dots correspond to models with f_b  = 0.5. The upper set of dots is for simulations of models that have e_b = 0.8, the middle set is for simulations of models that have e_b = 0.5, and the lower set is for simulations of models that have e_b = 0.2.

The quadrupole approximation made in deriving equation (5) is more accurate for larger planet orbital radii a. To test the accuracy of the prograde analytic solution, we consider simulations with a close-in planet. In Fig. 7, we plot equation (5) in the prograde case for e_b = 0.2 (blue line), 0.5 (black line), and 0.8 (red line) as a function of a. The points show simulation results for numerically determined stationary inclinations. The dot colours correspond to the binary eccentricity of the analytic lines. The orbit of the planet is unstable if a is less than about 2.3a_b. Thus, the innermost dots that we acquire by our simulations are for a planet at a = 2.3 (black line) and 2.4 (red and blue lines). The simulations are in very good agreement with the analytically results at large a, but deviate somewhat at smaller values of a. One indication of the deviation is that the time average eccentricity of the planet is about 0.1 for a = 2.5 a_b, rather than the value of 0 expected in the quadrupole approximation. The average eccentricity decreases with radius. At a = 3.0 a_b, it is about 0.04 and at a = 5.0 a_b, it is about 0.01.

Figure 7.

Comparison of the analytic solution given by equation (5) in the prograde case (solid lines) with simulation results (dots) for the stationary tilt i_s of the planet relative to the binary as a function of the semimajor axis of the planet a. The solid lines have e_b  = 0.2 (blue line), 0.5 (black line), and 0.8 (red line). The upper set of dots is for simulations of models that have e_b = 0.8, the middle set is for simulations of models that have e_b = 0.5, and the lower set is for simulations of models that have e_b = 0.2.

The black dotted lines in Figs 8 and 9 plot the prograde stationary inclinations i_s from our simulations and the magenta lines in the same figures plot the analytic solution for the prograde stationary inclination given by equation (5) as a function of binary eccentricity e_b for all parameters fixed except the binary eccentricity. Fig. 8 shows the results for a planet at |$r=5\, a_{\rm b}$| while Fig. 9 refers to the same system but with the planet at |$r=20\, a_{\rm b}$|⁠. The upper panels and bottom panels of two figures show the critical angles for an equal mass binary and a binary with f_b = 0.1, respectively. The left-hand and right-hand panels of the two figures show results for different planet masses. The black dotted lines are in very good agreement with the analytic results. The prograde stationary inclination values for the low-mass planet are rather insensitive to the location of the planet or f_b and so the lines look similar and are in the range of 80°–90°, since j is small (see Fig. 6). On the other hand, the stationary inclination is sensitive to e_b for the high-mass planet. The stationary inclination angle is smaller for smaller binary eccentricity and smaller binary mass fraction (larger j).

Figure 8.

Stationary inclination (i_s), critical minimum inclination (i_min) for libration, and critical maximal inclination for libration as a function of the binary eccentricity with the planet orbiting at r = 5a_b with different binary mass fractions f_b  = 0.5 (upper panels) and 0.1 (lower panels) for the lower mass planet (left-hand panels) and the high-mass planet (right-hand panels). The dotted lines plot the results of numerical simulations, while the solid lines are from the analytic model. The green dotted lines show the boundary between the prograde circulating and librating orbits while the blue dotted lines show the boundary between librating and retrograde circulating orbits. The black dotted lines show i_s obtained from our simulations. The magenta lines plot the analytic solutions for i_s from equation (5). The red lines plot the analytic solutions for i_min from equation (10) and the green solid lines plot the analytic solutions for i_min from equation (11).

Figure 9.

Same as Fig. 8 except that the planet is orbiting at r = 20a_b and the blue dotted line in the lower right panel shows the boundary between the polar libration and the crescent orbits region for the high-mass planet models.

3.1.2 Retrograde stationary inclination

The retrograde stationary inclination is given by equation (5), where we take the negative square root. The solid lines in Fig. 10 show the analytical solutions for e_b = 0.2 (blue line), 0.5 (black line), and 0.8 (red line) as a function of angular momentum ratio j. The retrograde stationary inclination, i_s, decreases monotonically with increasing j. However, the behaviour with binary eccentricity is more complicated, as we discuss further below. The six dots whose colours correspond to the e_b values for the analytic curves that represent i_s of Models G1–H3. The simulations are in very good agreement with the analytic predictions.

Figure 10.

Comparison of the retrograde analytic solution given by equation (5) with simulation results of Models G1–H3 for the retrograde stationary tilt i_s of the planet relative to the binary as a function of planet-to-binary angular momentum ratio j with binary eccentricity e_b = 0.2 (blue line), 0.5 (black line), and 0.8 (red line). The six dots represent the simulation results with j = 1 and j = 2. The curves reach i_s = 180° at j = j_cr given by equation (4).

There are two major qualitative differences between the prograde and retrograde stationary orbits. In the prograde case, for any value of planet-to-binary angular momentum ratio j there is a stationary inclination value about which there are nested librating orbits in the icos ϕ– isin ϕ plane (e.g. Figs 1 and 7). But in the retrograde case, we see from Fig. 10 that i_s reaches 180° for j = j_cr given by equation (4). For j < j_cr there are no stationary non-coplanar librating orbits, as discussed in Section 2.4.

The second difference between the prograde and retrograde stationary orbits is that for any fixed j, the stationary inclination angle i_s increases monotonically with e_b in the prograde case, but not generally in the retrograde case. The difference is seen by comparing Figs 7 (prograde) and 10 (retrograde) in that the curves in the prograde case for different e_b do not intersect, except at j = 0 and i_s = 90° that is the upper limit of prograde tilts, while in the retrograde case they do intersect and cross.

To confirm this crossing in the retrograde case, recognize that the intersection implies that i_s is independent of e_b at some fixed value of j. This condition can be written as

$$\begin{eqnarray*} \frac{\partial \cos {i_{\rm s}}}{\partial e_{\rm b}} \Big |_j=0, \end{eqnarray*}$$

(6)

where cos i_s is given by equation (5) in the retrograde case. This condition has an analytic solution for the intersection point at j = j_int and i_s = i_int, where

$$\begin{eqnarray*} j_{\rm int} &=& \frac{2}{\sqrt{3}}, \end{eqnarray*}$$

(7)

$$\begin{eqnarray*} i_{\rm int} &=& 150^\circ , \end{eqnarray*}$$

(8)

which is in agreement with the point of intersection in Fig. 10. For fixed j greater (smaller) than j_int, the stationary angle i_s increases with decreasing (increasing) binary eccentricity e_b.

3.2 Critical inclinations for libration

The critical minimum inclination angle between the prograde circulating and librating orbits can be determined analytically. We follow the description in section 3.4 of Martin & Lubow (2019). There are two branches, one with lower j and one with higher j, based on the sign of the parameter χ that is defined as

$$\begin{eqnarray*} \chi = e_{\rm b}^2-2(1-e_{\rm b}^2)j(2j+\cos i)\, . \end{eqnarray*}$$

(9)

The minimum tilt angle occurs where ϕ = 90°. If χ > 0, the minimum tilt angle for libration to occur is

$$\begin{eqnarray*} \cos i_{\rm min} = \frac{\sqrt{5}e_{\rm b}\sqrt{4e_{\rm b}^2-4j^2(1-e_{\rm b}^2)+1}-2j(1-e_{\rm b}^2)}{1+4e_{\rm b}^2}, \end{eqnarray*}$$

(10)

while if χ < 0, we have

$$\begin{eqnarray*} \cos i_{\rm min} = \frac{\sqrt{(1-e_{\rm b}^2)(1+4e_{\rm b}^2+60(1-e_{\rm b}^2)j^2)}-(1-e_{\rm b}^2)}{10(1-e_{\rm b}^2)j}\, . \end{eqnarray*}$$

(11)

The green and red solid lines in Figs 8 and 9 plot the analytic solutions of i_min obtained in equation (10) (green segments for χ < 0) and equation (11) (red segments for χ > 0) as a function of e_b. We also determined critical angles numerically in our simulations. The green dotted lines and blue dotted lines in Figs 8 and 9 show the simulation results for minimum and maximum inclinations, respectively. Figs 8 and 9 are determined for librating orbits of a low-mass planet (⁠|$m_{\rm p} = 0.001\, m_{\rm b}$|⁠) and a high-mass planet (⁠|$m_{\rm p} = 0.01\, m_{\rm b}$|⁠), respectively. Note that the dotted lines are shorter in the low binary mass fraction plots because some of the orbits are unstable in our simulations. However, the analytic solutions cover the entire range of binary eccentricities. The analytic solutions and the three-body simulations are in very good agreement. The comparison shows that both branches of the i_min analytic solutions (red and green) agree well with the simulations.

For a low-mass planet, the value for i_min is rather insensitive to the location of the planet or f_b and so the curves look similar in each panel. This insensitivity can be understood by the fact that such models have small values of angular momentum ratio j. In the limit that j goes to zero, we have from equation (9) that χ > 0 for e_b > 0. Therefore, i_min is given by equation (10) and the plotted curves should be nearly entirely red, rather than green. In this limit that j goes to zero, equation (10) reduces to

$$\begin{eqnarray*} \cos i_{\rm min}= \frac{\sqrt{5} e_{\rm b}}{\sqrt{1+4e_{\rm b}^2}}. \end{eqnarray*}$$

(12)

This equation implies that i_min decreases monotonically from 90° to 0° as e_b increases from 0 to 1, as we find in the low-mass planet plots.

In the limit that j is large, we have that χ < 0 in equation (9) and the entire curve for i_min should be green, as we find in the lower right panel of Fig. 9. In that limit, equation (11) applies and

$$\begin{eqnarray*} \cos i_{\rm min} = \sqrt{\frac{3}{5}}, \end{eqnarray*}$$

(13)

independent of e_b, as given by equation (40) in Martin & Lubow (2019). This minimum angle for libration is the same as the so-called Kozai–Lidov angle of ≃39.2° (Kozai 1962; Lidov 1962).

Therefore, in this high j limit, i_min should be constant and lie along the χ < 0 (green) branch, independent of e_b for e_b < 1. In the lower right panel of Fig. 9 (the highest j panel), we find that the green line is roughly what we predict. (There is a small red region close to e_b = 1 that is not visible in the plot.) The upper left panel in this figure has a smaller j value for a given e_b than in the lower right panel. In going from the former panel to the latter panel, we see that the behaviour of i_min is approaching the expectations of equation (13).

The lower right panel of Fig. 9 shows the results for the simulations with the same parameters as those in the upper right panel, but with a smaller binary mass fraction. In the smaller binary mass fraction case, there are no retrograde circulating orbits, only the crescent shaped orbits described in Fig. 4. Thus, the maximum libration angle (the blue dotted line) shows a different trend compared to the other parameters, which have retrograde precessing orbits.

4 DISCUSSION AND CONCLUSIONS

In this paper, we investigated the orbital evolution of a misaligned circular orbit planet with non-zero mass around an eccentric orbit binary by means of numerical simulations. The planet and binary interact gravitationally and the orbits of both vary in time. In particular, both undergo nodal precession in the inertial frame. In our suite of three-body simulations, we consider a low-mass planet with m_p = 0.001m_b at |$r=5\, a_{\rm b}$| and a high-mass planet with m_p = 0.01m_b at |$r=5\, a_{\rm b}$| and at |$r=20\, a_{\rm b}$| along with some even higher angular momentum third bodies. We considered different values of the eccentricity of the binary, e_b, its mass fraction, f_b, and the planet’s initial inclination i. To map out the possible orbits in these systems, we concentrated on numerically determining the transitions between the orbit families (circulating and librating for both prograde and retrograde orbits). In addition, we determined the stationary orbits for which the relative tilt and nodal phase between the planet orbit and binary orbit are constant in time.

For a very small planet mass, there are two stationary orbital states: coplanar and polar. In the polar state, the stationary planet-to-binary tilt is 90° and the angular momentum of the planet is along the binary eccentricity vector. The stationary states in the case that the planet mass is non-zero is a generalization to the polar state, but with the relative orientation not being perpendicular and the nodal phase not being constant in time in the inertial frame.

Equations (5), (9), (10), and (11) predict that the only parameters that control the planet-to-binary tilt for the stationary orbit and the minimum tilt for the transition from circulation to libration are the binary eccentricity e_b and the ratio of the planet-to-binary angular momentum j. Our numerical results agree with this prediction. Other parameters, such as the binary mass fraction f_b, only cause changes in these angles through their dependence on e_b or j. For example, in Fig. 6 we see that the stationary angle depends only on j for fixed e_b for different values of binary mass fraction f_b. In addition, the general agreement between the simulations and analytic predictions across a range of parameters implies that this dependence holds. These angles are also related to the evolution of discs. Simulations by Martin & Lubow (2019) suggest that a prograde disc approaches the stationary angle given by equation (5), if we consider j to represent the disc-to-binary angular momentum ratio.

The numerical results agree very well with the analytic equations for the stationary and minimum libration tilts given in Martin & Lubow (2019; see Figs 6, 8, 9, and 10). These analytic equations are based on the quadrupole approximation for the secular binary gravitational field (Farago & Laskar 2010). We find numerically that this approximation holds well, even for orbits that are fairly close to the binary ∼3a_b that is near the orbital radius where instability sets in (see Fig. 7).

As predicted analytically, the main effect of increasing angular momentum ratio j for fixed e_b is to monotonically decrease the relative planet-to-binary stationary tilt i_s in the prograde case (where i_s ≤ 90°) (see Fig. 6). In addition, this stationary tilt increases with increasing e_b for fixed j in the prograde case.

The behaviour of the stationary tilt in the retrograde case is more complicated, but agrees with the analytic predictions given in Martin & Lubow (2019). In this case, the stationary inclination for non-coplanar orbits decreases with increasing j for fixed e_b as in the prograde case. But i_s changes from increasing with e_b to decreasing at |$j=2/\sqrt{3}$|⁠. In addition, there are no non-coplanar stationary orbits below a certain j value, denoted by j_cr, that depends on e_b (see Fig. 10 and equation 4). This property does not hold in the prograde case.

Another difference between the prograde and retrograde cases is the topology of librating orbits. In the prograde case, the librating orbits are always nested in a icos ϕ–isin ϕ phase plane about point with i = i_s and ϕ = ±90° (e.g. red and cyan lines in Fig. 1). In the retrograde case, for j > j_cr described above, librating orbits are not fully nested within each other and do not orbit about the stationary point in the icos ϕ–isin ϕ phase plane (e.g. magenta lines in Fig. 5). In this phase plane, prograde librating orbits have an oval shape, while retrograde librating orbits have a crescent shape.

The variation of e_b in time is significantly larger for a higher mass planet. In the simulation with initial conditions e_b = 0.2, f_b = 0.1, and i = 170°, the binary eccentricity can be excited to values very close to 1 (see Fig. 1). This behaviour is similar to what occurs with Kozai–Lidov oscillations in which the outer object is a planet (e.g. Naoz 2016).

The recent release of data from Gaia allowed us to better characterize the kinematics of the eccentric equal mass close binary HD 106906 in the Lower Centaurus Crux group (Bailey et al. 2014). This system hosts both a wide asymmetric debris disc and a planetary-mass companion. The formation mechanism of the planet and the stability of its orbit are still under debate (Rodet et al. 2017; De Rosa & Kalas 2019). The binary components have masses |$M_1=1.37\, \rm M_\odot$| and |$M_2=1.34\, \rm M_\odot$|⁠. The orbital period is |$49.2\, \rm d$| and the binary eccentricity is e = 0.67 (De Rosa & Kalas 2019). The planet, HD 106906b, was directly imaged with a projected separation of 738 au and was found to be oriented at 21° from the position angle of the disc mid-plane, suggesting that the planet’s orbit is not coplanar with the system (Kalas et al. 2015). The planet mass was inferred to be |$m_{\rm p}=11\pm 2\, M_{\rm J}$| (Bailey et al. 2014). While the orbital properties of the planet are uncertain, if we assume that the semimajor axis is |$738\, \rm au$| and the eccentricity is zero, the ratio of the angular momentum of the planet to the angular momentum of the binary is j = 1.13. According to our analytical calculations, the stationary inclination is i_s = 73.5° and the minimum critical angle between the prograde circulating and the librating orbit regions is i_min = 41.6°.

The expected orbital properties of a high-inclination circumbinary planet depend on the mass of the gas disc in which it forms, and when the planet forms. Martin & Lubow (2019) showed that prograde discs of high mass align to the generalized polar state at tilts that are less than 90°. A giant planet that formed in such a disc would be expected to open a gap and not remain coplanar with the disc (Lubow & Martin 2016; Pierens & Nelson 2018). Such a massive planet would undergo libration oscillations in its orbit about an inclination of less than 90°, even after the disc has dispersed. As the disc loses mass, the inclination of the generalized polar state moves closer towards 90°. Thus, the time-scale of the disc dispersal may affect the final inclination of debris left over from the gas disc. A low-mass gas disc or massive disc that is dispersed on a sufficiently long time-scale forms a debris disc that orbits close to polar, as is the case for 99 Herculis, that is 3° away from polar alignment (Kennedy et al. 2012). A planet (e.g. an Earth-like planet) formed from the resulting polar debris disc would lie on a stationary polar orbit. A debris disc that is not close to polar or coplanar would be subject to violent collisions due to nodal differential precession.

ACKNOWLEDGEMENTS

Computer support was provided by UNLV’s National Supercomputing Center. CC acknowledges support from a UNLV graduate assistantship. We acknowledge support from NASA through grants NNX17AB96G and 80NSSC19K0443. Simulations in this paper made use of the rebound code which can be downloaded freely at http://github.com/hannorein/rebound.