Polar planets around highly eccentric binaries are the most stable

ABSTRACT

We study the orbital stability of a non-zero mass, close-in circular orbit planet around an eccentric orbit binary for various initial values of the binary eccentricity, binary mass fraction, planet mass, planet semimajor axis, and planet inclination by means of numerical simulations that cover 5 × 10⁴ binary orbits. For small binary eccentricity, the stable orbits that extend closest to the binary (most stable orbits) are nearly retrograde and circulating. For high binary eccentricity, the most stable orbits are highly inclined and librate near the so-called generalized polar orbit which is a stationary orbit that is fixed in the frame of the binary orbit. For more extreme mass ratio binaries, there is a greater variation in the size of the stability region (defined by initial orbital radius and inclination) with planet mass and initial inclination, especially for low binary eccentricity. For low binary eccentricity, inclined planet orbits may be unstable even at large orbital radii (separation |${\gt}5 \, a_{\rm b}$|⁠). The escape time for an unstable planet is generally shorter around an equal mass binary compared with an unequal mass binary. Our results have implications for circumbinary planet formation and evolution and will be helpful for understanding future circumbinary planet observations.

methods: analytical, methods: numerical, celestial mechanics, binaries: general, planetary systems

1 INTRODUCTION

Misaligned circumbinary discs are observed to be common in newborn binary systems (e.g. Chiang & Murray-Clay 2004; Winn et al. 2004; Capelo et al. 2012; Kennedy et al. 2012; Brinch et al. 2016; Kennedy et al. 2019). This may be a result of the chaotic accretion process during star formation (Bate et al. 2003; Bate 2018) or stellar flybys (Clarke & Pringle 1993; Cuello et al. 2019; Nealon, Cuello & Alexander 2020). Recent observations show that there is a break between alignment and misalignment of circumbinary discs with the central binary at orbital periods of about |$30 \, \rm d$|⁠. About 68 per cent of short period binaries (period |${\lt}20\, \rm d$|⁠) have aligned discs (within 3°), while those with longer orbital period show a larger range of inclinations and binary eccentricities (Czekala et al. 2019). Since circumbinary planets form within a circumbinary disc, we expect their initial alignments to be similar to the disc.

Misaligned discs should be more common around eccentric orbit binaries (Martin & Lubow 2017). Eccentric orbit binaries occur at longer binary periods. Main-sequence binaries are observed to be on circular orbits for periods up to about 8 d due to tidal circularization (e.g. Raghavan et al. 2010). Binary eccentricities appear somewhat limited for binaries with periods of less than about 30 d, likely due to the effects of stellar tidal dissipation. Of the 10 circumbinary planets observed through transits by the Kepler mission, all of them are in orbits that are close to coplanar to their host binary orbit (Orosz et al. 2012a,b; Welsh et al. 2012; Kostov et al. 2014; Welsh et al. 2015; Kostov et al. 2016; Li, Holman & Tao 2016). All the Kepler transit detected circumbinary planets orbit around a low-eccentricity binary, except for Kepler-34b with a binary eccentricity of 0.52 (Welsh et al. 2012; Kley & Haghighipour 2015). However, the observed coplanarity may be a selection effect because these planets are found around binaries with smaller orbital periods and generally low binary eccentricity. In addition, the Kepler data may contain planets transiting non-eclipsing binaries and more inclined circumbinary planets may be found in the future (Martin & Triaud 2014).

Binaries with longer orbital periods are expected to host planets with a wide range of inclinations (relative to the binary orbital plane) because their planet-forming discs have a wide range of inclinations. But planets around longer period binaries are harder to detect by transit methods because the transit probability is smaller and because the planet orbital period is longer. Eclipse timing variations of the binary may be able to distinguish polar planets (those that are perpendicular to the binary orbital plane) from coplanar planets (Zhang & Fabrycky 2019). Furthermore, there is evidence based on binary eclipse timing variations for a highly misaligned circumbinary planet around the binary star KIC 5095269 (Getley et al. 2017). KIC 5095269b has an estimated mass of |$7.70 \, M_{\rm J}$| and orbits a binary with a moderate eccentricity of 0.26 (Getley et al. 2017). However, results from Borkovits et al. (2016) suggest that the binary eccentricity is lower e_b = 0.05 with a planet orbit inclination of 40°.

For a test (massless) particle that represents a low-mass planet in orbit around a circular orbit binary, its angular momentum vector precesses around the angular momentum vector of the binary with constant tilt. Since the nodal precession rate does not change sign, the longitude of the ascending node fully circulates over 360°. The particle’s orbit is circulating and may be either prograde or retrograde relative to the binary, depending upon the initial particle inclination. On the other hand, for a binary with non-zero eccentricity, a circumbinary test particle orbit may undergo libration if the test particle has a sufficiently large initial inclination. The longitude of the ascending node covers a limited range of angles, less than 360°, and the nodal precession rate changes sign for a librating orbit. In this case, the angular momentum vector of the planet precesses about the binary eccentricity vector (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 a planet to undergo libration decreases with increasing binary eccentricity and so a planet with even a small inclination may undergo libration in a highly eccentric binary. We define the stationary orbit of a test particle to be the orbit which does not change secularly in time in the inertial frame. That is, its orbit averaged tilt and longitude of ascending node are constant over long time-scales. A particle orbit that is coplanar with a binary has this property. In addition, there is a stationary inclination state in which the particle angular momentum vector lies along the binary semimajor axis (eccentricity vector). In this so-called polar aligned state, the orbital plane of the particle is perpendicular to the binary orbital plane.

For a particle with non-zero mass, the binary orbit is no longer fixed. In this case, the stationary orbit, called the generalized polar orbit, refers to an orbit that is fixed in a frame comoves with the binary orbit. That is, the orbit averaged tilt and longitude of ascending node of the generalized stationary orbit do not vary on long (secular) time-scales in a frame defined by the orbit-averaged directions of the binary eccentricity and angular momentum. As the mass of the planet (particle) increases, the orbital inclination for the generalized polar aligned state decreases (Chen et al. 2019; Martin & Lubow 2019).

A circumbinary disc may also undergo either circulation or libration depending on its initial inclination, binary eccentricity, and the disc angular momentum. For a disc that is sufficiently warm and radially limited, it may undergo solid body precession (Larwood & Papaloizou 1997). Dissipation within the disc leads to alignment, either coplanar alignment to the binary orbital plane (Papaloizou & Terquem 1995; Lubow & Ogilvie 2000; Nixon et al. 2011; Nixon, King & Price 2012; Facchini, Lodato & Price 2013; Foucart & Lai 2013, 2014; Lodato & Facchini 2013) or polar alignment (Aly et al. 2015; Martin & Lubow 2017, 2018; Zanazzi & Lai 2018; Franchini, Lubow & Martin 2019a; Smallwood & et al. 2019). For an eccentric orbit binary, a disc undergoes tilt oscillations as it aligns, even if it is aligning to coplanar (Smallwood et al. 2019). However, if the lifetime of the disc is shorter than the disc alignment time-scale of a planet-forming disc, then planets may form in a misaligned disc that is neither coplanar nor polar. Thus, a giant planet may have a misaligned orbit around a binary system. The planet and the disc may not remain coplanar with each other once the planet opens a gap in the disc (Picogna & Marzari 2015; Lubow & Martin 2016; Martin et al. 2016; Franchini, Martin & Lubow 2019b).

In this work, we study the orbital stability of misaligned close-in planets around eccentric binaries for different binary and planet parameters. Holman & Wiegert (1999), Sutherland & Fabrycky (2016), and Quarles et al. (2018) considered the orbital evolution and stability of a test particle on an initially coplanar circular orbit about an eccentric binary. Doolin & Blundell (2011) extended these test particle results to the non-coplanar case (see also Quarles & Lissauer 2016; Hong et al. 2019). In Chen et al. (2019), we further extended these studies to investigate the evolution of an inclined circumbinary planet with mass on an initially circular orbit. The orbital behaviour of the planet is in good agreement with the analytic results of Martin & Lubow (2019) that is based on the quadrupole approximation for the binary potential (Farago & Laskar 2010). The agreement breaks down for close-in planets where the quadrupole approximation is less accurate. In the octupole approximation, these orbits have only been studied for test particles and not planets with mass (Naoz et al. 2017; Vinson & Chiang 2018).

We found that there is a generalized polar orbit whose inclination angle depends on the planet mass. This orbit is stationary in a frame that precesses with the binary. Recently, Cuello & Giuppone (2019) and Giuppone & Cuello (2019) included an analysis of the orbital dynamics of massive P-type planets around eccentric binaries. They considered up to Jupiter mass planets with an orbital inclination of 90°. This configuration does not provide a generalized polar (stationary) orbit, as is consistent with their findings. We describe the initial conditions for the three-body simulations and show the orbital stability maps in Section 2. In Section 3, we discuss the implications of our results and draw our conclusions.

2 THREE-BODY SIMULATIONS

In this section, we first describe the simulation set-up and the parameter space that we explore. We then show the stability maps that describe the stable circumbinary orbits. Finally, we consider the time-scale on which the unstable orbits become unstable.

2.1 Simulation set-up and parameter space explored

To study the stability of a third body orbiting around an eccentric binary star system, we use a whfast integrator that is a second-order symplectic Wisdom Holman integrator with 11th-order symplectic correctors in the N-body simulation package, rebound (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 for which the central binary has components of mass m₁ and m₂ with total mass m_b = m₁ + m₂. The mass fraction of the binary is defined with f_b = m₂/m_b. The separation of the binary is a_b, the eccentricity of the binary is e_b, and the orbital period of the binary is T_b.

The planet is initially in a Keplerian orbit around the centre of mass of the binary and has mass m_p. Its orbit is defined by six orbital elements: the semimajor axis a, inclination i relative to the binary orbital plane, eccentricity e, longitude of the ascending node ϕ measured from the binary semimajor axis, argument of periapsis ω, and true anomaly ν. The orbit of the planet is initially circular and so initially e = 0 and ω = 0. We choose ϕ = 90° or ϕ = 0° initially in our suites of simulations. In Table 1, we describe the sampling of the phase space for how we vary a, i, and ν. Note that the binary orbit is not fixed since the binary feels the gravity of the massive third body.

Table 1.

The initial conditions for the planet orbits in the three-body simulations. The planet is initially in a circular Keplerian orbit about the centre of mass of the binary. The initial planet eccentricity is e = 0, the argument of periapsis is ω = 0, and the longitude of the ascending node is ϕ = 90°. Column 1 shows the name of the orbital parameter. Column 2 shows the minimum value in the phase space. Column 3 shows the maximum value in the phase space. Column 4 shows the spacing in the phase space.

Orbital element	Minimum value	Maximum value	Δ
a	1.5 ab	6 ab	0.05 ab
i	0	π	π/80
ν	0	5π/3	π/3

In this work, we are interested in the stability of circumbinary planet orbits that are close to the binary. We consider orbits for which the initial semimajor axis 1.5a_b ≤ a ≤ 6a_b over times of 5 × 10⁴T_b and over the full range of initial inclinations. We define the orbit as unstable once at least one of three criteria is met. Instability occurs first, if the eccentricity of the planet becomes large, e > 1.0, so that the planet is not bound to the binary; second, if the semimajor axis of the planet increases significantly, |$a \gt 10\, a_{\rm b}$|⁠; or third, if the semimajor axis of the planet becomes very small, a < a_b (see e.g. Quarles et al. 2019). We define the escape time to be the time at which one of these criteria is satisfied. If the planet does not meet any of the three criteria for instability within the time of 5 × 10⁴T_b, then we classify it as stable.

In order to classify the orbital motion of the planet, we work in a frame relative to the instantaneous eccentricity and angular momentum vectors of the binary (⁠|$\boldsymbol {e_{\rm b}}$| and |$\boldsymbol {l_{\rm b}}$|⁠). The frame has the three axes along unit vectors |$\boldsymbol {\hat{e}_{\rm b}}$|⁠, |$\boldsymbol {\hat{l}_{\rm b}}\times \boldsymbol {\hat{e}_{\rm b}}$|⁠, and |$\boldsymbol {\hat{l}_{\rm b}}$|⁠. For planet angular momentum |$\boldsymbol {l_{\rm p}}$|⁠, the inclination of the planet’s orbital plane relative to the binary orbital plane is given by

$$\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 planet. The inclination of the binary relative to the total angular momentum |$\boldsymbol {l}$| is

$$\begin{eqnarray} \beta_{\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, the phase angle (longitude of ascending node) of the particle in the same frame of reference is given by

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

(3)

(Note that this equation corrects equation 3 in Chen et al. 2019 that contains some typos.)

There are four types of orbits in our stability maps. These are shown in the icos ϕ – isin ϕ phase space in fig. 1 in Chen et al. (2019). There are two types of circulating orbits, meaning the planet angular momentum precesses about the binary angular momentum vector. The green lines represent prograde orbits that display retrograde (clockwise) precession, dϕ/dt < 0, and the blue lines represent retrograde orbits that display prograde (counterclockwise) precession, dϕ/dt > 0. The librating orbits, meaning those where the planet angular momentum vector precesses about the direction of angular momentum of the stationary orbit, are shown in the cyan and red lines. The stationary inclination for a test particle is 90°, but the generalized stationary state for a planet with mass is at lower inclination. The inclination at the centre of the librating orbits is the stationary inclination, i_s. The red orbits have initial inclination i < i_s, while the cyan orbits have initially i > i_s. These librating orbits display prograde precession in ϕ. The centre is at i = i_s ≈ 90° for a low-mass planet case and it is <90° for a high-mass planet case (see equation 5 in Chen et al. 2019).

Figure 1.

Orbital stability as a function of initial planet separation a and initial inclination i for initial ϕ = 90°. The binary is equal mass, f_b = 0.5, and the initial binary eccentricity is e_b = 0.2 (first row), 0.5 (second row), and 0.8 (third row). The third body has mass |$m_{\rm p} = 10^{-3}\, m_{\rm b}$| (first column), |$5\times 10^{-3}\, m_{\rm b}$| (second column), and |$10^{-2}\, m_{\rm b}$| (third column). The green and blue pixels represent prograde and retrograde circulating orbits, respectively. The red pixels represent librating orbits with initial inclination i < i_s while the purple pixels represent librating orbits with initial inclination i > i_s. Each pixel represents simulations for six different values of the true anomaly and the darker the colour, the more stable orbits. The white pixels represent unstable orbits.

Table 2 lists the values of the initial binary eccentricity, binary mass fraction, and planet mass for all of the simulations that we consider. We explore three binary eccentricities, e_b = 0.2, 0.5, and 0.8, two binary mass fractions f_b = 0.5 and 0.1, and three different planet masses, |$m_{\rm p} = 0.001\, m_{\rm b}$|⁠, |$0.005\, m_{\rm b}$|⁠, and |$0.01\, m_{\rm b}$|⁠.

Table 2.

Parameters of the simulations. The first column contains the name of the Model, the second and third columns indicate the binary eccentricity and mass fraction. The fourth column represents the mass of the planet in units of m_b.

Model	eb	fb	mp (mb)
A1	0.2	0.5	0.001
A2	0.2	0.1	0.001
A3	0.2	0.5	0.005
A4	0.2	0.1	0.005
A5	0.2	0.5	0.01
A6	0.2	0.1	0.01
B1	0.5	0.5	0.001
B2	0.5	0.1	0.001
B3	0.5	0.5	0.005
B4	0.5	0.1	0.005
B5	0.5	0.5	0.01
B6	0.5	0.1	0.01
C1	0.8	0.5	0.001
C2	0.8	0.1	0.001
C3	0.8	0.5	0.005
C4	0.8	0.1	0.005
C5	0.8	0.5	0.01
C6	0.8	0.1	0.01

In the case with initial ϕ = 90° there are always librating orbits. This can be seen for example from fig. 1 in Chen et al. (2019). For test particles all possible circulating and librating orbits with fixed binary eccentricity are covered by sampling along the vertical line in the icos ϕ − isin ϕ, corresponding to ϕ = 90°. In the non-zero mass planet case, the vertical line does not sample all possible orbits at fixed initial binary eccentricity because the binary eccentricity varies in time. For example, for the case with initial ϕ = 90°, two planet orbits with different initial inclinations but the same initial binary eccentricity have different binary eccentricity by the time they reach ϕ = 0°. We are considering stability maps where the initial binary eccentricity is the same across all orbits so it is necessary to consider orbits at different initial phase angles. For each model, we consider two initial values for the nodal phase angle, ϕ = 90° and ϕ = 0°. The case of initial ϕ = 90° is most favourable for producing librating orbits, while the case of initial ϕ = 0° does not produce any librating orbits.

2.2 Stability maps for initial ϕ = 90°

We first consider the stability maps for equal mass binaries and then for a binary mass fraction of f_b = 0.1.

2.2.1 Equal mass binary

Fig. 1 shows stability maps for equal mass binary models with initial binary eccentricity e_b = 0.2 (models A1, A3, and A5, first row), models with binary eccentricity e_b = 0.5 (models B1, B3, and B5, second row), and models with binary eccentricity e_b = 0.8 (models C1, C3, and C5, third row). In each panel, we vary the initial separation, a, and the initial inclination, i. For each combination of these two parameters considered, we vary the true anomaly ν from 0° to 300° with an interval of Δ = 60°. Thus, each pixel in Fig. 1 represents six simulations in total.

The colour of each pixel represents the type of planet orbit. We follow the colour representation in Chen et al. (2019), except we use purple instead of cyan. Thus, the prograde circulating orbits are green and the retrograde circulating orbits are blue. The librating orbits with initial inclination i < i_s are red and those with initial inclination i > i_s are purple. Notice that one pixel can only display one colour for the six simulations in each pixel. In some cases, there may be two types of orbits in the same pixel. In such cases, the colour refers to the type of orbit that has the largest number of the six stable orbits. The darker the colour, the larger the number of stable orbits. White pixels indicate unstable orbits. A lighter non-white pixel indicates that some contributing orbits are unstable.

As seen in Fig. 1, the prograde circulating (green) orbits are the generally least stable orbits for the equal mass binary, meaning that the minimum initial semimajor axis a for stable orbits is generally larger than for the other orbit types. The critical inclination below which the orbit is prograde circulating (green) (i.e. the inclination at red–green boundary) is higher for smaller binary eccentricity and higher planet mass. The critical inclination does not vary much with planet separation for the small separations considered, as expected in the analytic model of Martin & Lubow (2019). For these green orbits, the minimum stable initial planet semimajor axis does not change much with planet inclination or planet mass, but is more strongly affected by the binary eccentricity. The higher the binary eccentricity, the more unstable the prograde circulating (green) orbits are.

In Fig. 1, the retrograde circulating (blue) orbits are the most stable orbits for low binary eccentricity e ≲ 0.5. The critical inclination above which the orbit is retrograde circulating (blue) (i.e. inclination at the purple–blue boundary) is lower for smaller binary eccentricity and lower planet mass. These circulating retrograde orbits are generally stable closer to the binary than the circulating prograde orbits. But there is significant dependence of stability on the initial planet orbit inclination with orbits starting closer i = 180° being more stable. For the equal mass binary case shown here, there is not much difference in stability of the retrograde circulating orbits across different values of the planet mass (i.e. across a row in Fig. 1).

The librating orbits are the most stable orbits for high binary eccentricity. The stationary inclination decreases with increasing planet mass. That is, the inclination i_s at the red–purple boundary decreases with increasing planet mass, as expected in the analytic model of Martin & Lubow (2019). The range of inclinations for which there are librating orbits with initial i < i_s (red orbits) decreases with planet mass while the range of inclinations for which librating orbits with initial i > i_s (purple orbits) increases with planet mass. The most stable orbits at high eccentricity occur near the stationary inclination. Generally, the larger the difference in the initial planet inclination from the stationary inclination |i − i_s|, the more unstable the librating planet orbits.

Comparing with fig. 14 in Doolin & Blundell (2011), we find that our results are in general agreement with their test particle simulations, since the mass of the planet has not affected our results much for an equal mass binary. However, the stationary inclination, i_s, is always 90° in their simulations because the test particle does not have angular momentum to exchange with the binary system. In addition, we have considered higher binary eccentricity than their simulations that cover up to e_b = 0.6. Consequently, we see more clearly that the near polar librating orbits are the most stable at high binary eccentricity for an equal mass binary.

2.2.2 Binary with mass fraction f_b = 0.1

Fig. 2 shows the stability maps for binary mass fraction f_b = 0.1 for binary eccentricity e_b = 0.2 (models A2, A4, and A6, first row), binary eccentricity e_b = 0.5 (models B2, B4, and B6, second row), and binary eccentricity e_b = 0.8 (models C2, C4, and C6, third row). The range of angles for which the planet is librating is larger compared with the equal mass binary case, however there is much more instability for low binary mass fraction. The effect of increasing the mass of the planet is much more significant for f_b = 0.1 than it is for an equal mass binary, as a consequence of the planet mass being closer to the binary secondary mass.

Figure 2.

Same as Fig. 1 except the binary mass fraction is f_b = 0.1.

The prograde circulating (green) orbits are the least stable orbits for high binary eccentricity, e_b ≥ 0.5, similar to the equal mass binary case. The innermost stable prograde orbits have |$a \simeq 2.3\, a_{\rm b}$| for the simulation with e_b = 0.2. This radius does not change much with the mass of planet for low binary eccentricity. However, the innermost stable orbits are much farther out for high binary eccentricity, the separation is |$a \simeq 3.5 \, a_{\rm b}$| for binary eccentricity e_b = 0.8. The critical inclination below which the orbit is prograde and circulating (i.e. inclination at the green–red boundary) increases by about 3° as the planet mass increases for e_b = 0.2, while for e_b = 0.8, it increases by about 15° with planet mass increasing from |$0.001\, m_{\rm b}$| up to |$0.01\, m_{\rm p}$|⁠. These results are consistent with analytic solutions of the green–red lines in the bottom panels of fig. 8 in Chen et al. (2019) that show that a low-mass planet has higher i_crit for small e_b and lower i_crit for large e_b. The boundary between prograde circulating and the librating orbits is unstable even if the planet is at large separation, |$a \gt 5 \, a_{\rm b}$|⁠, in the high-mass model (m_p = 0.01 m_b) with low binary eccentricity e_b = 0.2 (Model A6). However, with increasing binary eccentricity, Model B6 and Model C6 show that the boundary becomes more stable.

For low binary eccentricity e_b ≤ 0.5, the retrograde circulating (blue) orbits are again the most stable orbits. The innermost stable retrograde region is at separation |$a \simeq 1.5 \, a_{\rm b}$| for e_b = 0.2 and |$a \simeq 3.2 \, a_{\rm b}$| for e_b = 0.8. The critical inclination above which the orbit is retrograde circulating (i.e. inclination at the purple–blue boundary) is lower for smaller e_b and lower m_p. It increases by about 20° with increasing m_p for all the binary eccentricities considered here. For the high e_b model, the critical inclination of the retrograde circulating orbits increases very close to 180° for |$m_{\rm p} = 0.01\, m_{\rm b}$|⁠. This is consistent with blue dotted lines in the two bottom panels of fig. 8 in Chen et al. (2019) that show the critical angles of the low-mass planet and the high-mass planet increase by about 20° as the binary eccentricity increases from e_b = 0.2 to 0.8 in the binary system with f_b = 0.1.

As in the equal mass binary case, the librating orbits are the most stable orbits for high binary eccentricity. But in this case of an unequal mass binary, the difference in stability between librating and circulating orbits is greater than in the equal mass binary case, especially at high binary eccentricity. The innermost stable libration region is at |$a \simeq 2.5 \, a_{\rm b}$| for e_b = 0.2 and |$a \simeq 2.3 \, a_{\rm b}$| for e_b = 0.8 for the models with m_p = 0.01 m_b. The stationary inclination i_s decreases with increasing m_p. Thus, the inclination of the transition from red to purple orbits i_s decreases with increasing planet mass. For the small e_b models (e_b = 0.2), the purple libration region expands with increasing planet mass m_p, while the red libration region becomes very small. Thus, there are no red librating orbits in the icos ϕ–isin ϕ phase plot of Model D1 in fig. 3 in Chen et al. (2019).

For the test particle orbits considered in Doolin & Blundell (2011), the librating orbits in the stability maps appear close to symmetric about the stationary inclination i = i_s = 90°. However, the symmetry is broken for a planet with mass. The asymmetry is larger for lower e_b and higher m_p. For the low m_p models, even for a planet at large radius (a > 5a_b), the purple librating orbits with i = 100° < i < 120° are unstable while the red libration region is stable. This is consistent with the icos ϕ–isin ϕ phase plot of Model B1 in fig. 1 in Chen et al. (2019) where there were no stable orbits for initial i in this region. On the other hand, for the high m_p model at small binary eccentricity e_b = 0.2, the orbit of the planet at large radius (a > 5a_b) with i = 60° < i < 80° is unstable. This is consistent with the icos ϕ–isin ϕ phase plot of Model D1 in fig. 3 in Chen et al. (2019) where there were no red stable librating orbits for initial i in this range.

Overall, for initial ϕ = 90°, at low binary eccentricity, the polar orbits are generally more unstable for the low-mass fraction binary compared to the case of an equal mass binary. The retrograde circulating orbits are the most stable orbits for low e_b and the innermost region of the stable retrograde orbit can extend down to |$a = 1.5 \, a_{\rm b}$|⁠. The mass of the planet does not affect the stability map much for equal mass binary, but it has a significant impact for the unequal mass binary, particularly for low binary eccentricity. For high binary eccentricity, the librating orbits are generally the most stable, particularly for inclinations close to the generalized polar angle i_s.

Cuello & Giuppone (2019) reported dynamical maps for polar orbits for massive planets with chaotic regions. They also find that polar circumbinary planets around eccentric binaries are less stable for low binary mass ratios (see their fig. 8). Giuppone & Cuello (2019) investigated the orbital stability of a planet with i = 90° orbiting a binary with e_b = 0.5 for various values of f_b and m_p (see e.g. their fig. 2). Their results for m_p = 0.001 |$\rm M_{\odot }$| show that the planet is less stable for low binary mass ratios and they are consistent with our simulation results in the middle left panel of Fig. 2.

As a check on our results, we applied the mean exponential growth of nearby orbits (MEGNO) indicator value 〈Y〉 that can identify whether the orbit is chaotic or quasi-stable (Cincotta & Simó 2000). If 〈Y〉 converges to a value 2 over time, the orbit is stable, while if 〈Y〉 diverges linearly in time, the orbit is chaotic. To test how the MEGNO indicator compares with our results, we ran the most extremely eccentric and high-mass case (Model C6) shown in the lower right of Fig. 2 using the high accuracy IAS15 integrator. The resulting stability map in Fig. 3 shows most of the stable orbits in the lower right of Fig. 2 have 〈Y〉 ≃ 2. Therefore, the MEGNO results are quite similar and support our conclusion that the generalized polar orbits are the most stable for high binary eccentricity.

2.3 Stability maps for initial ϕ = 0°

In the case with initial ϕ = 0° there are no librating orbits. This can be seen for example from fig. 1 in Chen et al. (2019). Along the horizontal line in the icos ϕ − isin ϕ phase plot there are no librating orbits.

Fig. 4 shows the stability maps for initial ϕ = 0° and an equal mass binary. Unlike Fig. 1 for the ϕ = 90° case, we see that the librating region is absent. The circulating regions are enlarged with similar features. Initial inclination of 90° means that the initial angular momentum vector of the planet orbit is perpendicular to both the binary angular momentum vector and the binary eccentricity vector, in the direction of the |$\boldsymbol {j}\times \boldsymbol {e_{\rm b}}$| vector.

Figure 3.

Same as Model C6 in Fig. 2 with the high accuracy IAS15 integrator and based on the MEGNO indicator. The colour points are for MEGNO indicator 〈Y〉 between 1.5 and 2.4 at a time of 5 × 10⁴T_b that suggests orbital stability.

Figure 4.

Same as Fig. 1 except initial ϕ = 0°.

Fig. 5 shows the stability maps for initial ϕ = 0° and binary mass fraction f_b = 0.1. The orbits are generally less stable compared to the equal mass binary case. There is a more sensitive dependence on the mass of the planet in this case, as we found for the initial ϕ = 90° case. The stability of the orbits in the ϕ = 0° case is quite similar to the stability for the circulating orbits in the ϕ = 90° case, but extends to a wider range of initial inclinations.

Figure 5.

Same as Fig. 2 except initial ϕ = 0°.

The critical initial inclination above which the orbits are circulating and retrograde (inclination at the green–blue boundary) increases with increasing planet mass. For the equal mass binary simulations, the critical inclination of the retrograde orbits with a high-mass planet is about 100°. For small binary mass fraction, the critical initial inclination is larger than in the corresponding equal mass binary case and increases with increasing binary eccentricity. The critical initial inclination with high planet mass and low binary fraction increases to about 140° for e_b = 0.8.

2.4 Circumbinary planet escape time

Figs 6 and 7 show the escape times of unstable orbits, t_esc. The panels in Figs 6 and 7 correspond to the panels in Figs 1 and 2, respectively. Each pixel in the density plot represents the average t_esc of an orbit. Darker pixels correspond to shorter escape times and a white pixel represents a stable orbit.

Figure 6.

Average escape time, t_esc, as a function of planet semimajor axis a for a binary system with binary mass fraction f_b = 0.5 and initial binary eccentricity e_b = 0.2 (first row), 0.5 (second row), and 0.8 (third row). The third body has mass |$m_{\rm p} = 10^{-3}\, m_{\rm b}$| (first column), |$5\times 10^{-3}\, m_{\rm b}$| (second column), and |$10^{-2}\, m_{\rm b}$| (third column).

Figure 7.

Same as Fig. 6 except the binary mass fraction is f_b = 0.1.

We find that a system with smaller f_b has more long-lived unstable orbits (more light pixels) and a system with equal mass fraction models has more short-lived unstable orbits (more dark pixels), similar to that described in Doolin & Blundell (2011) for the test particle orbits.

Moreover, orbits with lower initial inclination are generally more short-lived and the number of short-lived unstable orbits with high initial inclination increases with increasing binary eccentricity. Orbits with initial inclination close to the stationary inclination are more long-lived orbits even in the innermost region especially for models with the small f_b. For more extreme mass ratio binaries, the binary torques due to resonances are weaker and therefore the unstable planet orbits evolve more slowly.

Finally, we do not observe a significant difference for unstable orbits between systems with different planet mass in the equal mass binary models. On the other hand, there are more long-lived unstable orbits with increasing planet mass in the lower f_b models especially for e_b = 0.2.

3 DISCUSSION AND CONCLUSIONS

In this paper, we have investigated the orbital stability of a misaligned initially circular orbit close-in planet with non-zero mass around an eccentric orbit binary by means of numerical simulations. In our suite of simulations, we have considered planets with masses m_p = 0.001, 0.005, and |$0.01\, m_{\rm b}$| around a binary with mass fraction f_b = 0.1 and 0.5 and we sample different orbits by varying the initial semimajor axis, a, inclination, i, and true anomaly, ν. We consider two values for the initial nodal phase angle, ϕ = 0° and ϕ = 90°.

In general, we find that circumbinary planet orbits are stable over times of 5 × 10⁴T_b for initial orbital radii a that are greater than about 2.2a_b to about 4a_b, although there are cases where the instability extends to larger radii. The values depend on the planet and binary parameters.

The results show that at high binary eccentricity, planet orbits near the generalized polar orbits are the most stable type of planet orbit. In particular, they are more stable than prograde and retrograde coplanar orbits. The enhanced stability of these generalized polar orbits over coplanar orbits is most apparent for more a extreme mass ratio, high eccentricity binary, as seen in the bottom right of Fig. 2.

The range of radii covered by stable nearly polar orbits increases with binary eccentricity. For high binary eccentricity, the resonant binary Lindblad torque on a low-mass polar gas disc decreases with increasing binary eccentricity (e.g. Miranda & Lai 2015; Lubow & Martin 2018; Franchini et al. 2019a). As a result, the disc inner edge of a polar disc extends closer to the binary centre of mass than it does for a prograde coplanar disc. Analogous results appear in the planet orbit case in that the radii for stable orbits extend closer to the binary centre of mass for polar orbits than is the case for coplanar orbits. For the polar case involving a high eccentricity binary, the potential due to the binary is at each instant of time is nearly axisymmetric in the plane of the planet orbit. Consequently, the torque exerted on the planet perpendicular to its orbit plane is relatively weak compared to a lower eccentricity case. Just the opposite behaviour would be expected for a coplanar planet. That is, with increasing binary eccentricity the binary potential becomes more non-axisymmetric in the plane of the planet orbit leading to a larger torque. Such considerations support the simulation results we find that the polar orbits are more stable than coplanar ones.

The orbital alignment of circumbinary planets depends upon the final alignment of the circumbinary disc. A low-mass disc may evolve towards either coplanar alignment or polar alignment at i = 90°. However, high-mass discs may align to the generalized polar state at tilts that are less than 90° (Martin & Lubow 2019). Assuming the disc evolves through instantaneously stationary (generalized polar) configurations, the disc may move towards a tilt of 90° as it loses mass. Therefore, the time-scale of the disc dispersal plays a role in the final inclination of a debris disc or protoplanet. Disc–planet interactions must also be taken into account. A giant planet that opens a gap may not remain coplanar with a misaligned disc (Lubow & Martin 2016; Pierens & Nelson 2018; Franchini et al. 2019b). Thus, a massive planet could undergo libration oscillations in its orbit even after the disc has dispersed.

Our results have implications on the possible orbits in which circumbinary planets can reside. The results in Figs 1 and 2 suggest that for highly eccentric orbit binaries, planets found inward of about 3.5a_b can reside on only misaligned librating orbits including polar orbits. Planets found inward of about 2.5a_b around low-eccentricity binaries can reside on only circulating retrograde orbits. Planets at larger separations |${\ga}4 a_{\rm b}$| can be stable for both librating and circulating orbits over a wide range of planet orbit inclinations. That is, orbital stability considerations alone do not provide a constraint on the possible planet orbits, if the planet is far enough from the binary. But based on disc simulations, we expect there to be a preference for highly inclined and polar planets around high-eccentricity binaries (e.g. Martin & Lubow 2018). Whether a disc evolves to polar depends on its initial inclination and nodal phase. The detailed occurrence frequency of polar disc configurations and therefore polar planets then depends on the initial disc misalignment distribution that is not known. But binaries with inclined circumbinary discs are found to occur for binary orbital periods longer than about 30 d (Czekala et al. 2019).

A promising technique for finding polar orbit planets is to use binary eclipse timing variations (ETVs). To date, planet detections using ETVs based on Kepler data have been reported for the inclined planet KIC 5095269b (Borkovits et al. 2016; Getley et al. 2017). Two additional systems KIC 07177553 and KIC 7821010 may also have an inclined planetary mass object orbiting around a binary (Borkovits et al. 2016). Since the ETV technique is becoming well developed, it is likely that additional misaligned and polar circumbinary planets will be found by the future observations with TESS or PLATO (Zhang & Fabrycky 2019). We believe that more inclined planets will be found, especially around longer period, eccentric orbit binaries.

ACKNOWLEDGEMENTS

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