© ESO 2019

1 Introduction

Planets appear as common by-products of the star formation process within molecular clouds (Chiang & Youdin 2010). At early evolutionary stages (class 0/I), the fraction of stellar binaries and higher-order multiple systems is remarkably high: between 30 and 70% (Connelley et al. 2008; Chen et al. 2013). Consequently, protoplanetary discs can be severely affected by binary formation and stellar flybys (e.g. Bate 2018; Pfalzner 2013). Given the stochasticity of these processes, highly asymmetrical discs naturally form on timescales of the order of a few hundred kyr. Therefore, planet formation is also expected to occur in a broad variety of misaligned discs (Bate 2018; Cuello et al. 2019a). Currently, about a hundred planets have been detected in multiple-star systems (Martin 2018). There is, however, a striking tension with the more than 3700 planets detected¹ around single stars (Batalha et al. 2013).

Planets in binary systems can either be circumbinary (P-type, also known as Tatooines) or circumstellar (S-type). The seeming lack of P-type planets has been ascribed to observational biases and dynamical processes. On the one hand, these systems are intrinsically difficult to observe through radial velocities (Eggenberger & Udry 2007; Wright et al. 2012) and transit methods (Martin & Triaud 2014). For transits in non-coplanar systems, the arbitrary orientation of the planet chord coupled to the fact that most binaries do not eclipse renders P-type planet detection extremely challenging. On the other hand, recent stellar-tidal evolution models of short-period binaries show that the binary orbital period increases with time (Fleming et al. 2018). This translates into a larger region of dynamical instability around the binary, which could explain the lower frequency of P-type planets compared to S-type planets. Based on the available data, typical planets found around binary stellar systems have a radius of the order of 10 Earth radii and orbital periods of about 160 days (Martin 2018). For instance, around a binary system with total mass equal to 1 M_⊙, this corresponds to an orbit with a semi-major axis of approximately 0.35 au.

Previous theoretical studies suggested that the disc-binary interaction, namely the gravitational torque, should align the disc with the binary orbital plane on timescales shorter than the disc lifespan (Foucart & Lai 2013). Therefore, until recently, misaligned circumbinary discs (CBDs) were considered exotic. However, the observations of highly non-coplanar systems such as 99 Herculis (Kennedy et al. 2012), IRS 43 (Brinch et al. 2016), GG Tau (Cazzoletti et al. 2017), and HD 142527 (Avenhaus et al. 2017) suggest otherwise. For example, in the multiple system of GG Tau, Aly et al. (2018) recently explored various initial conditions to explain the disc misalignment with respect to the binary. Even more recently, Kennedy et al. (2019) reported the very first observation of a circumbinary gas-rich disc in a polar configuration. These discoveries are rather puzzling and deserve to be explained.

In the field of secular dynamics, the seminal work by Ziglin (1975) presented an integrable approximation for the orbits of planets around binary systems. Later, based on numerical calculations of circumbinary polar particles, Verrier & Evans (2009) reported the coupling between the inclination and the node. Then, Farago & Laskar (2010) studied the circumbinary elliptical restricted three-body problem, giving an averaged quadrupolar Hamiltonian. In this way the authors showed that equilibrium configurations can be reached at high inclinations in the restricted and massive three-body problem. Then, relying on this analysis, it is possible to define the separatrix around the equilibrium points.

A more thorough analysis of the stability of massless particles was done by Doolin & Blundell (2011) where they consider different binary eccentricities and mass ratios between the components. This approach allows us to unambiguously identify prograde, retrograde, and polar orbits. In addition, the authors study the validity of the quadrupolar Hamiltonian approximation when the distance to the binary centre of mass is not sufficiently small. Furthermore, based on the study of hierarchical triple systems, several other studies investigated the circumbinary polar orbits for the restricted problem considering the octupole Hamiltonian (Ford et al. 2000; Naoz et al. 2013; Li et al. 2014). Interestingly, for near polar configurations, the inner binary can significantly excite the orbital eccentricity (up to about 0.3 in some cases). This excitation is weakly dependent on the distance to the binary. More recently, Naoz et al. (2017) and Zanardi et al. (2018) studied the evolution of polar orbits considering additional effects such as general relativity terms. Last but not least, the effect of the binary on the CBD can have dramatic effects on the disc alignment as described below.

Given the recent works on polar alignment and CBD dynamics (Sect. 2), it is interesting to explore further the conditions under which polar alignment of misaligned CBDs occurs (Sect. 3). Then, assuming planet formation happens in these discs, we characterise the orbits of single P-type planets around eccentric and inclined binaries (Sect. 4). This allows us to constrain the binary eccentricity and mass ratio of systems around which we expect to observe polar circumbinary planets in the future. We discuss our findings in Sect. 5 and finally conclude in Sect. 6.

2 Summary of previous works

In the following we list the series of works that have previously addressed the topic of misaligned CBDs. We pay particular attention to nodal libration, polar CBDs, and alignment conditions.

2.1 Nodal libration, binary parameters, and stability

Doolin & Blundell (2011) numerically investigated the nodal libration mechanism in the longitude of the ascending node and in the inclination for P-type particles around eccentric binaries. This phenomenon occurs for highly inclined orbits and for a wide range of binary mass ratios. Polar orbits were found to be remarkably stable at distances of two to three times the binary semi-major axis, which is not the case for coplanar configurations. Moreover, between the islands of stability, the authors report vertical striations of instability and hypothesise that they might be due to orbital resonances (see their Fig. 14, and Sect. 4.1). More recently, Zanazzi & Lai (2018) derived a simple analytical criterion based on the binary eccentricity and initial disc orientation that assesses whether the inclination of the CBD is expected to librate or not. It is worth noting that this criterion is symmetric with respect to the disc orientation: it does not distinguish prograde from retrograde configurations. However, as we show in Sect. 3.3, this condition is not sufficient in case of symmetry breaking.

2.2 Polar alignment of CBDs

In Aly et al. (2015) and Martin & Lubow (2017) it is demonstrated through numerical simulations how inclined gaseous CBDs around eccentric binaries can tidally evolve towards polar configuration with respect to the binary orbital plane. In this case, the dissipative torque acting on the disc is key to damping the nodal libration. This mechanism is investigated further in Zanazzi & Lai (2018) and Lubow & Martin (2018) where an analytical framework is provided to describe the polar alignment of CBDs, taking into account the gaseous viscosity of the disc.

2.3 Flybys and generation of polar alignment conditions

Stellar flybys, also known as tidal encounters, provide a natural way to tilt and twist a CBD. More specifically, a disc initially aligned with the binary orbital plane can be tilted by several tens of degrees and twisted by roughly 90° after a flybyof a massive star (Cuello et al. 2019a, Fig. D6). In particular, the smaller the pericentre distance and the more massive the perturber, the higher the resulting disc misalignment. It is worth noting that this flyby must happen in the early stages of disc formation such that enough time should remain for the disc to reach the polar configuration. In this scenario, inclined andretrograde encounters are the most efficient for disc warping (Xiang-Gruess 2016). Alternatively, misaligned accretion onto the binary, as in Bate (2018), provides another mechanism to reach appropriate conditions for polar alignment.

2.4 CBDparameters and alignment timescales

More recently, Martin & Lubow (2018, hereafter ML18) explored the alignment of CBDs in a broad variety of configurations through hydrodynamical simulations. Specifically, they studied the dependence of the disc viscosity, temperature, size, binary mass ratio, orbital eccentricity, and inclination on the disc evolution. The authors find that polar alignment occurs for a wide range of binary–disc parameters, especially for low mass discs. Their results broadly agree with the expectation of linear theory. Interestingly, for large binary separations and/or for low mass ratios, the alignment timescale might be longer than the expected CBD lifetime. If this is the case, the CBD might remain on a misaligned configuration, neither (anti-)aligned nor polar. This suggests that highly inclined discs and planets may exist around eccentric binaries (Zanazzi & Lai 2018).

3 Circumbinary disc dynamics and polar alignment

We first summarise the conditions under which we expect polar alignment. Then we present our SPH simulations, and finally describe the unexpected symmetry breaking observed.

3.1 Binary parameters and polar orbits

We consider a binary system with total mass M and individual masses M₁ and M₂. We call q = M₂∕M₁ the binary mass ratio. To describe the motion of a P-type particle around the binary we use the following Jacobi orbital elements: semi-major axis a, eccentricity e, inclination i (with respect to the binary orbital plane), mean longitude λ, longitude of pericentre ϖ (alternatively argument of pericentre ω), and longitude of the ascending node Ω. The sub-index B is used when referring to the binary orbit. Without loss of generality, we set initial conditions {λ_B = 0°, ϖ_B = 0°, Ω_B = 0°}. The angles of a given particle are measured from the direction of the binary pericentre. Finally, we set our fiducial system M₁ = M₂ = 0.5 M_⊙ (i.e. q = 1), a_B = 0.1 au, and e_B = 0.5. This allows us to directly compare our results with those in Martin & Lubow (2017). Regarding the dynamics of P-type particles, the phase space depends on Ω with two equilibrium points at {Ω = ±90°, i = 90°}, where both angles librate (see Doolin & Blundell 2011, for a detailed analysis). Interestingly, for any given couple of (Ω, i), it is possible to assert whether i librates around ± 90° or not by calculating the separatrix F (Farago & Laskar 2010): $$F=\sqrt{\frac{1-e_{\text{B}}^2}{1-5e_{\text{B}}^2{\cos }^2\Omega +4e_B^2}}.$$(1)

Specifically, its orbit is called polar if a particle fulfils the following condition: $$\arcsin F<i<\pi -\arcsin F.$$(2)

In Fig. 1 we show the separatrix in the plane (Ω, i) for the fiducial parameters considered. The initial conditions inside the separatrix limits correspond to the region of polar orbits for this given configuration. We recall, as described in Sect. 2.2, that the higher the value of e_B, the larger the region.

To sum up, in order to be polar, a P-type particle needs to have a longitude of the ascending node close to 90° and be orbiting around an eccentric binary on a highly inclined orbit (see Fig. 1). When these conditions are met and in the absence of gas, then the particle’s inclination and the ascending node librate around 90°.

3.2 Misaligned circumbinary disc simulations

In this section, we assume that the chaotic and violent stellar cluster environment has produced an initially tilted CBD (see Sect. 2.3). We computed its resulting alignment, considering different initial disc inclinations (i₀) and twists (Ω₀). For simplicity, we only considered the set of fiducial parameters for the binary.

We ran 3D hydrodynamical simulations of CBDs by means of the PHANTOM Smoothed Particle Hydrodynamics (SPH) code (Price et al. 2018). The SPH method (see e.g. Price 2012) is well suited for misaligned CBDs simulations because there is no preferred geometry and angular momentum is conserved to the accuracy of the time-stepping scheme, irrespective of the plane of the binary orbit. Both stars are represented by sink particles with accretion radii equal to 0.25 a_B, which interact with the gas via gravity and accretion (Bate et al. 1995). We used 10⁶ gas particles to model the disc, which is a resolution high enough to avoid artificially high numerical viscosities². The inner and outer radii of the disc were set to R_in = 2 a_B and R_out = 5 a_B, respectively.The surface density was initially set to a power law profile equal to Σ ∝ R^−3∕2. The disc total mass is equal to 0.001 M_⊙. Given this low value, we neglected the disc self-gravity in our calculations. Furthermore, we assumed that the disc is locally isothermal (i.e. that the sound speed follows c_s ∝ r^−3∕4 and H∕R = 0.1 at R = R_in). For further details on the PHANTOM disc set-up (see Price et al. 2018). Finally, we adopted a mean Shakura–Sunyaev disc viscosity α_SS ≈ 0.01.

Our set of PHANTOM simulations is given in Table 1. We only change the disc inclination (i₀) and its argument of ascendingnode (Ω₀). For completeness, these disc initial conditions are plotted as numbered dots in Fig. 1. Based on the linear theory in Zanazzi & Lai (2018), the disc is expected to become polar in 2 and 4, align in 1 and 6, and anti-align in 3 and 5, with respect to the binary plane. From the practical point of view, we decompose the disc in 500 annuli and compute the tilt and the twist for each one of them. There is little difference between the disc warping in the inner and in the outer regions, except for r3 (see Sect. 3.3). Therefore, unless stated otherwise, we compute the tilt (i) and the twist (Ω) at r = 3 a_B from the PHANTOM outputs.

In Fig. 2 we plot the evolution of the tilt over time for 1000 binary orbital periods. We obtain the expected alignment for all the runs but r3. On the one hand, in r1 and r6 the disc inclination remains prograde and is slowly damped, while Ω circulates (not shown). In r5 the disc anti-aligns with respect to the binary (i = 180°), while Ω circulates. On the other hand, in r2, r3, and r4 the disc becomes polar (i = 90°) and Ω librates (instead of circulating). Remarkably, r2 and r4 are in phase opposition and the oscillations are quickly damped over time (compared to r1 and r6). Also, for low inclinations (r5 and r6) the disc becomes coplanar more quickly in the retrograde configuration (r5) compared to the prograde (r6) because the inner cavity, opened by resonant torques, is smaller for retrograde configurations (Miranda & Lai 2015; Nixon & Lubow 2015). Therefore, the disc responds more quickly to the binary torques. The anomalous behaviour of r3 is detailed below.

Fig. 1 Separatrix in the (i, Ω) plane constructed from Eq. (1) for eB = 0.5. Initial conditions inside the separatrix limits correspond to polar orbits 2 and 4, whereas those outside correspond to prograde (i < 90°: 1, 6) and retrograde (i > 90°: 3, 5). The initial conditions of the SPH simulations of Table 1 are identified with numbers.

Fig. 2 Tilt evolution over time, i(t), in terms of binary orbital periods. Runs: r1 (solid blue), r2 (dashed blue), r3 (solid red), r4 (dashed red), r5 (dotted red), and r6 (dotted blue). For comparison with r3, we also plot the symmetric curve of r1 with respect to i  =   90° (solid green). r3 shows an anomalous behaviour because the disc breaks in two discs after ~ 150 orbits, which then evolve towards polar alignment (see Sect. 3.3).

3.3 Unexpected symmetry breaking in r3

In Fig. 3 we show the CBD in r3 for three evolutionary stages: 0, 165, and 1000 binary orbital periods. Given the initial parameters of r3, the disc should anti-align with respect to the binary (see Fig. 1). However, after half a tilt oscillation (~ 150 binary periods), the disc suddenly fulfils the condition for polar alignment and its tilt starts to oscillate around 90°. Additionally, its twist starts to librate instead of circulating. Both oscillations are damped, as expected for the viscous evolution of a gaseous disc (ML18). Hence, this shows that the symmetry between prograde and retrograde orbits is broken in some specific cases (e.g. r3).

In Fig. 4, we show the evolution of the ratio between the outer (a = 5 a_B) and the inner (a = 2 a_B) regions of the CBD in r3. This quantity indicates whether the disc is behaving as a solid body (~ 1) or not. By construction, this ratio is initially equal to 1. However, after a hundred binary orbits, it starts to oscillate between 0.8 and 1.2 due to the disc breaking, which is caused by the binary torque (see middle row in Fig. 3). We observe that the inner disc regions undergo a fast polar alignment, which eventually causes the entire CBD to align at i = 90°. After 400 binary orbits the entire disc has reached polar alignment and behaves as a solid body again (i.e. the tilt ratio is close to 1). This phenomenon is representative of the symmetry breaking between prograde and retrograde configurations.

The polar alignment in r3 is anomalous in the sense that it cannot be predicted by the analytical theory. More specifically, the expression of Λ in Eq. (7) in Zanazzi & Lai (2018), which shows whether the CBDs is expected to become polar or not, is symmetric for prograde and retrograde orbits (see e.g. their Fig. 2). Therefore, fulfilling the polar alignment criterion at the very beginning is not a necessary condition for polar alignment. The cavity size and further non-linear effects play a key role in determining the final CBD alignment. Finally, as in ML18, we find that during the simulations the variations in the binary orbital parameters are only of the order of a few percent.

Polar CBDs might then be more common than previously thought. Therefore, since planet formation is expected to occur within these misaligned discs, it is relevant to characterise the planetary orbits around eccentric and inclined binaries.

4 P-type planet orbital characterisation

We focus on the stability and evolution of single circumbinary P-type planets (a > a_B) for a wide variety of binaries. To do so, we construct a regular 2D mesh of initial conditions with a pair of orbital parameters (e.g. i and e), while the remaining four parameters (e.g. a, λ, ω, and Ω) are initially set at nominal values. Then, we solve the three-body equations of motion numerically by means of a double-precision Bulirsch-Stoer integrator with tolerance 10⁻¹².

4.1 Dynamics of an inclined massless planet

First, we focus on the case of massless particles. A thorough study of the survival times for a regular combination of binary eccentricity 0 < e_B < 0.9 and mass ratio 0.2 < q < 1 in the (a, i) plane can be found in Doolin & Blundell (2011). Here, we extend their analysis by focusing on the features observed in the dynamical maps in order to better characterise inclined orbits. In particular, we use further structure indicators.

In Fig. 5, we show the dynamical maps in the (a, i) plane for Ω =0° (leftmost panels) and Ω = 90° (rightmost panels). Each orbit is integrated for 80 000 binary periods. If the particle either collides with one of the stars or escapes from the system, then it is coloured in white. In the top row, the colour scale is proportional to the Mean Exponential Growth of Nearby Orbits (MEGNO) value ⟨Y ⟩. This indicator is particularly appropriate because it identifies the chaotic orbits at a low CPU-cost (Cincotta & Simó 2000). Regular orbits have ⟨Y ⟩ ~ 2 (blue), while chaotic and potentially unstable orbits have ⟨Y ⟩ > 2 (red/brown). In the bottom row, the colour scales correspond to the amplitude variation of i.

On the one hand, when Ω = 0° (first column in Fig. 5), the prograde orbits are stable only beyond 3.5 a_B, while retrograde orbits are stable closer to the binary system (a ~ 2.2 a_B). This occurs because the overlap of nearby resonances for retrograde orbits is of higher order compared to the case of prograde orbits (see Morais & Giuppone 2012). The vertical red stripes are associated with the 1 : N mean motion resonances (MMR), where eccentricity is efficiently excited. The orbits of these region have ⟨Y ⟩ ≫ 2, which indicates their chaotic nature. For prograde cases the 1 : 7 MMR delimits the stable orbits from the unstable ones at roughly a ~ 3.5 a_B. Alternatively, for retrograde cases this happens for the 1 : 4 MMR at a ~ 2.5 a_B. Interestingly, for high initial inclinations the orbits are stable for both prograde and retrograde configurations, even fairly close to 90°. It is worth noting that this does not occur for S-type configurations (see white gaps in Fig. 5 in Giuppone & Correia 2017).

On the other hand, when Ω = 90° (rightmost columns in Fig. 6), the stability maps change significantly and different stable regions appear. Specifically, there are stable orbits with high inclinations (45° ≤ i ≤ 135°) as close as retrograde orbits. Polar orbits are remarkably stable at small semi-major axis (a ~ 2.5 a_B). Furthermore, we analyse the impact of the binary mass ratio q on the stability of the orbits. The two rightmost columns of Fig. 5 show the stability maps for q = 0.5 and q = 0.2. As the mass ratio decreases, the MMRs become more important and severely deplete the polar regions (i ~ 90°). These vertical stripes in the e_B –q parameter space were first reported by Doolin & Blundell (2011). Here, we demonstrate that these peculiar patterns areindeed located at MMRs. In addition for non-equal-mass binaries the libration amplitude of i (and consequently of the nodes) increases for polar orbits (see Fig. 5, bottom panel). This effect correlates well with the chaotic nature of these orbits.

To further characterise the orbits, we construct additional dynamical maps in a regular (e_B, i) grid. In Fig. 6, we show these maps for a fixed semi-major axis equal to a = 4.5 a_B and Ω =90°. This corresponds to a regular region of motion according to Fig. 5. Again, this is done for different binary mass ratios: q = 1, 0.5, and 0.2. To visually identify the region of polar orbits, we overplot a continuous red line in Fig. 6 defined by Eq. (2). The regions with the highest eccentricity excitation (Δe≳ 0.5) or with e_B ≥ 0.945 are unstable (coloured in white). When q = 1 (leftmost panel) the dynamical maps exhibit low eccentricity amplitudes with MEGNO values indicating regular orbits. For lower values of q, horizontal stripes of unstable and chaotic orbits appear at i = i_N and i = i_c. The values of i_N ~ 73°and 134° correspond to the very narrow octupole resonance wherein ϖ + ϖ_B − 2Ω librates about 0°, while i_c ~ 46°and 107°, where the frequency $\dot{\varpi }~0$ at ω± 180° (see e.g. Gallardo et al. 2012; Vinson & Chiang 2018, for a detailed analysis of these configurations). As a rule of thumb, the lower the mass ratio q, the larger the size of the chaotic region at i = 90°.

We also observe a diamond-shaped feature centred at i = 90° and delimited by the separatrix with a vertex at e_B ~ 0.45. We checked that the regions where Δe ≳ 0.4 have high values of MEGNO (not shown) and are hence chaotic. In the bottom panel of Fig. 6, we show longitudinal cuts of the dynamical maps at i = 90°. The curves show the variations in eccentricity (Δe) and inclination (Δi) for different values of e_B. In the leftmost panel (q = 1) there is no discontinuity observed for Δe or Δi. However, for q = 0.5, we observe an abrupt change in Δe and Δi at e_B ~ 0.42. For q = 0.2, this occurs when 0.4 ≲ e_B ≲ 0.5. These curves clearly show that Δe and Δi are strongly correlated.

Fig. 3 Circumbinary disc in r3: the binary plane lies in the xy-plane and the disc is initially inclined (i0 = 120°, Ω0 = 0°). At 165 binary orbits, we observe that the tilt is different between the inner and the outer regions of the disc (cf. Fig. 4). This is interpreted as the disc being broken in two discs by the binary. The inner regions experience polar alignment, causing the whole disc to become polar after roughly 400 binary orbits (cf. Fig. 2). At 1000 binary orbits, the disc is almost perfectly polar.

Fig. 4 Tilt ratio between a = 5 aB and a = 2 aB in r3. When the disc behaves as a solid body, this ratio is ~1. The oscillations of this ratio (between 0.8 and 1.2) appear when the disc breaks in two misaligned discs. This occurs between approximately 150 and 400 binary orbits (see middle row in Fig. 3).

4.2 Orbital characterisation for a massless polar planet

In the following, we exclusively focus on polar orbits (i.e. those that are initially at Ω = 90°, i = 90°). Figure 7 shows the polar orbits in the (e_B, q) plane. The colour scale corresponds to Δe. High values of Δe ≳ 0.3 result into chaotic behaviour. We see that some of the orbits (white dots) are unstable in less than 10⁵ binary periods. This dynamical map allows us to conclude that binaries with 0.35 ≲ e_B ≲ 0.5 and 0.05 ≲ q ≲ 0.4 are un- likely to harbour polar planets. In other words, this corresponds to the forbidden region for polar type-P circumbinary planets.

4.3 Orbital characterisation for a massive polar planet

Here we set the eccentricity of the binary to e_B = 0.5 and vary themass ratio q, and the mass of the planet M_p for polar orbits. By doing so we aim to explore whether the stability depends on the planetary mass. In Fig. 8, we show the dynamical maps obtained for ΔΩ and for Δ(Ω − ϖ_B). For q > 0.7 and M_p< 10⁻⁵ M_⊙, the node librates (ΔΩ ~ 0°). For the same mass ratios but higher planetary masses we see that the node librates with higher amplitudes (ΔΩ > 90°). Interestingly, for massive planets M_p > 10⁻⁵ M_⊙ (i.e. giants) the orbit precesses regardless of the value of q. This occurs because the massive planet causes the binary to start precessing. In addition, we note that the value of Δ(Ω − ϖ_B) does not depend on the planetary mass. Hence, in this configuration the indicator ΔΩ allows us to discern whether the binary ellipse is fixed in time or not. Figure 8 (lower panel) shows that, as expected, both the orbits of a single giant planet and of an almost massless planet remain polar around the osculating binary ellipse (i.e. independently of the binary precession).

Finally, in Fig. 9 we show the evolution of four selected orbits identified with numbers in Fig. 8. Orbit (1) corresponds to an orbit of an almost massless (10⁻¹⁵ M_⊙, roughly the mass of a planetary embryo) around a binary with q = 0.73. Orbit (1) is to be compared with orbit (2), which corresponds to an embryo as well, but orbiting around a binary with q = 0.67. Orbit (2) exhibits secular excitation in e and Ω, and it is stable for at least 10⁸ binary orbits. Orbit (3) corresponds to a planet of mass 10⁻³ M_⊙ (~ Jupiter mass) orbiting around a binary with q = 0.67. We observe that the node circulates and that the secular excitations in Δe and Δi are of the same order as in orbit (2). Finally, orbit (4) corresponds to the orbit of a small moon (10⁻¹¹ M_⊙) around a binary with q = 0.2. In this case the excitation in e reaches values of 0.4, which may eventually trigger instability. In fact, according to the MEGNO indicator, orbits 1, 2, and 3 are regular (⟨Y ⟩~ 2), while orbit 4 is chaotic (⟨Y ⟩ > 2). Based on the varied behaviour of these characteristic orbits, we conclude that planets with high masses (M_p ≳ 1 × 10⁻⁵ M_⊙) cause the binary pericentre to precess, and that the most regular orbits are obtained around binaries with q > 0.7.

In summary, the most favourable conditions to obtain stable polar orbits for circumbinary planets are reached, at moderate to high inclinations, around equal-mass components and eccentric binaries for which Ω ~±90°. In addition, planets with masses above 10⁻⁵ M_⊙ strongly affect the binary orbit. Assuming that such a massive planet forms in the CBD, its gravitational perturbations are likely to render the surrounding disc material unstable. The study of this effect is beyond the scope of this work and will be more thoughtfullyexplored in the future.

5 Discussion

As mentioned in Sect. 2.3, the specific initial conditions for polar alignment can either be reached after a stellar flyby (Cuello et al. 2019a) or through misaligned accretion within the molecular cloud (Bate 2018). Therefore, it is reasonable to consider polar alignment of the CBD as a possible (and potentially likely) dynamical outcome (Aly et al. 2015; Zanazzi & Lai 2018; Lubow & Martin 2018). The presence of the gaseous disc around the binary is crucial in order to viscously damp the nodal libration mechanism of the CBD (see Sect 2.2.). This is indeed a necessary condition for polar alignment, which should be distinguished from polar oscillations.

The final disc configuration, however, is highly sensitive to the binary parameters. For example, increasing the binary eccentricity e_B (for a fixed Ω) widens the inclination range for which polar alignment is expected. Additionally, increasing the binary mass ratio q decreases the oscillation period (cf. Eq. (16) in Lubow & Martin 2018) and makes the disc behave less rigidly. This is precisely when symmetry breaking between prograde and retrograde configurations is expected to occur. Then, for high values of q and e_B, the CBD is more susceptible to breaking up with at least one of the two discs eventually becoming polar.

The disc parameters are equally important as shown by ML2018: a higher disc viscosity leads to a faster damping of the nodal libration mechanism. Also, cool discs could be more easily broken by the binary torque. When the disc precession timescale is shorter than the sound crossing time, then the disc breaks at roughly 3 a_B. Finally, the larger the disc, the longer the timescale for the oscillations, even though these oscillations are more rapidly damped.

In this context, the very recent detection of a gas-rich CBD in polar configuration by Kennedy et al. (2019) is of particular interest. In other cases, the presence of a highly inclined disc – perhaps in a polar configuration – has been proposed to explain asymmetric illumination patterns (Marino et al. 2015; Min et al. 2017; Facchini et al. 2018). In some cases, these shadows are able to trigger spiral formation in the outer disc (Montesinos et al. 2016; Montesinos & Cuello 2018), which are able to efficiently trap dust (Cuello et al. 2019b). Therefore, a shadow-casting inner disc is expected to affect the process of planet formation in the outer disc. Interestingly, the presence of shadows and spirals in an outer CBD could provide an indirect way to infer the presence of an inner CBD in a polar configuration.

Then, provided planet formation is efficient in these systems, polar and misaligned circumbinary P-type planets should be common around mildly eccentric (e ≤ 0.4) equal-mass binaries (q > 0.7). We recall that planet formation on polar orbits is more favourable when δe and δi are small since these conditions guarantee a small velocity dispersion. This kind of planet is stable at distances of a few times the binary semi-major axis (a ≳ 3a_B). For Ω ~ 90°, polar stability islands appear at even closer distances from the binary. These regions should be considered as sweet spots for polar Tatooine formation and evolution.

Interestingly, planets with masses above 10⁻⁵ M_⊙ significantly affect the binary orbit and exhibit fast node circulation (alternatively binary pericentre precession). Observational campaigns such as BEBOP (Martin et al. 2019) provide, and will continue to provide, useful constraints on their potential existence and frequency.

Last but not least, we also expect gaps at 1 : N resonances with the binary once the gas dissipates in the remaining circumbinary debris disc. In the absence of any viscous force, these gaps naturally appear for high values of ⟨Y ⟩ in Fig. 5. This effect has profound implications for planetesimal growth and planet evolution in more evolved CBDs, which will be investigated in a future work.

Fig. 5 Dynamical maps of inclined systems with M = 1 M⊙ (from left to right): {q = 1, Ω = 0°}, {q = 1, Ω = 90°}, {q = 0.5, Ω = 90°}, and {q = 0.2, Ω = 90°}. The top labels correspond to the period in terms of binary periods. The regions between the horizontal red lines correspond to real polar orbits (i.e. librating around i = 90°). Long-term integrations show that particles with initial conditions in the red/brown regions do not survive more than 105 periods. When Ω = 90°, stability islands appear at polar inclinations closer to the binary (between 2 and 3 aB). In the bottom panel, the colour scale shows the amplitude of libration of i. The continuous red lines are defined by Eq. (2) and correspond to the analytic limits between polar and prograde or retrograde orbits.

Fig. 6 Top panels: dynamical maps of inclined configurations with Ω = 90° with q = 1, 0.5, 0.2 (from left to right). The colour scale corresponds to the eccentricity variation Δe of the orbit. The continuous red lines are defined by Eq. (1) and correspond to the analytic limit between polar and prograde or retrograde orbits. Bottom panels: longitudinal cuts of the dynamical maps at i = 90°. The amplitude of Δe (Δi) is shown in red (blue). The grey shaded areas correspond to chaotic orbits with ⟨Y ⟩ ≫ 2.

Fig. 7 Dynamical map for polar orbits for massless planets with initial a = 4.5 aB in the (eB, q) plane. The colour scale corresponds to Δe. The orbits identified with white dots are unstable. Chaotic orbits identified with the MEGNO indicator are indicated with grey squares.

Fig. 8 Dynamical maps of polar orbits (i = 90°, Ω = 90°) for massive planets in the (q, Mp) plane with eB = 0.5. The colour scale in the upper and lower panels corresponds to ΔΩ and Δ(Ω − ϖB), respectively. Initially, the planet is located at the equilibrium point i = 90°, Ω = 90°. The four initial conditions identified with white circles correspond to the orbits shown in Fig. 9. Each system was integrated for 20 000 yr (~ 650 000 binary periods).

Fig. 9 Evolution of the individual polar orbits of Fig. 8: (1) blue, (2) green, (3) brown, and (4) yellow. The secular coupling between eccentricity, inclination, and node is easily identified. The periods are of roughly 5 × 104 binary periods, depending on the couple of values (MP, q). The evolution of giant planets (MP > 10−4, i.e. brown region in Fig. 9) exhibits the circulation of node. According to the MEGNO indicator, orbits 1, 2, and 3 are regular (⟨Y ⟩~ 2); instead, 4 is chaotic after ~100 000 orbital periods of the binary (⟨Y ⟩ > 2).

6 Conclusions

We first studied the evolution of misaligned CBDs through hydrodynamical simulations (see Sect. 3). Then, by means of N-body integrations,we thoroughly characterised the dynamics and stability of single P-type planets around eccentric binaries (see Sect. 4). The main results of this study can be summarised as follows:

• 1.

  The viscous torque exerted by the binary on the misaligned CBD is stronger for retrograde configurations compared to progradeones. Hence, for the same relative inclination with respect to i = 90°, a retrograde CBD is more likely to become polar compared to a prograde CBD. This kind of symmetry breaking can happen in the inner disc regions under specific initial conditions (e.g. our run 3). Due to tilt and twist oscillations, the polar alignment condition of the CBD (initially not fulfilled) can be reached after a few hundred orbits. When this happens, the CBD stops behaving as a solid body (see Figs. 3 and 4) and the disc breaks. This phenomenon is not captured by the linear theory and is more likely to occur around eccentric equal-mass binaries.

• 2.

  Our N-body simulations show that although “exotic” polar circumbinary P-type planets are stable on long timescales, namely over 10⁷ binary orbital periods. Based on stability criteria, we expect to find this kind of planet around binaries with mild eccentricities (e_B < 0.4) and for any value of the mass ratio q (see Fig. 7).

• 3.

  We find that above a certain planetary mass threshold (around 10⁻⁵ M_⊙ in Fig. 8), the massive polar Tatooine strongly affects the binary orbit, which then starts precessing. The impact of the binary precession and the gravitational perturbations of such a giant planet on a hypothetical disc of planetesimals will be the subject of a future work.

In conclusion, planet dynamics around misaligned and polar CBDs is dramatically different compared to similar coplanar configurations. We note that a lack of detections of polar circumbinary planets, mainly because of observational biases, does not constitute proof of non-existence. Future radial-velocity and transit surveys around binaries will help to constrain the fraction of polar, inclined, and retrograde P-type planets. Equal-mass binaries with eccentricities ~ 0.4 should be considered as the most promising targets to discover polar Tatooines. Thus, since HD 98800BaBe has a high eccentricity (e_B = 0.785 ± 0.005, Kennedy et al. 2019), this system does not constitute a good candidate for harbouring such a planet.

Acknowledgements

We thank the anonymous referee and the editor for very valuable comments and suggestions that have significantly improved our work. We also thank Cristian Beaugé, Agnese Da Rocha Rolo, and Antoine Rocher for useful discussions during this project. N.C. acknowledges financial support provided by FONDECYT grant 3 170 680. The authors acknowledge support from CONICYT project Basal AFB-170 002. The Geryon cluster at the Centro de Astro-Ingenieria UC was extensively used for the SPH calculations performed in this paper. BASAL CATA PFB-06, the Anillo ACT-86, FONDEQUIP AIC-57, and QUIMAL 130 008 provided funding for several improvements to the Geryon cluster. N-body computations were performed at Mulatona Cluster from CCAD-UNC, which is part of SNCAD-MinCyT, Argentina. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 823823. Figure 3 was made with SPLASH (Price 2007).

References

All Tables

All Figures

	Fig. 1 Separatrix in the (i, Ω) plane constructed from Eq. (1) for eB = 0.5. Initial conditions inside the separatrix limits correspond to polar orbits 2 and 4, whereas those outside correspond to prograde (i < 90°: 1, 6) and retrograde (i > 90°: 3, 5). The initial conditions of the SPH simulations of Table 1 are identified with numbers.
In the text

Fig. 2 Tilt evolution over time, i(t), in terms of binary orbital periods. Runs: r1 (solid blue), r2 (dashed blue), r3 (solid red), r4 (dashed red), r5 (dotted red), and r6 (dotted blue). For comparison with r3, we also plot the symmetric curve of r1 with respect to i  =   90° (solid green). r3 shows an anomalous behaviour because the disc breaks in two discs after ~ 150 orbits, which then evolve towards polar alignment (see Sect. 3.3).

In the text

Fig. 3 Circumbinary disc in r3: the binary plane lies in the xy-plane and the disc is initially inclined (i0 = 120°, Ω0 = 0°). At 165 binary orbits, we observe that the tilt is different between the inner and the outer regions of the disc (cf. Fig. 4). This is interpreted as the disc being broken in two discs by the binary. The inner regions experience polar alignment, causing the whole disc to become polar after roughly 400 binary orbits (cf. Fig. 2). At 1000 binary orbits, the disc is almost perfectly polar.

In the text

	Fig. 4 Tilt ratio between a = 5 aB and a = 2 aB in r3. When the disc behaves as a solid body, this ratio is ~1. The oscillations of this ratio (between 0.8 and 1.2) appear when the disc breaks in two misaligned discs. This occurs between approximately 150 and 400 binary orbits (see middle row in Fig. 3).
In the text

Fig. 5 Dynamical maps of inclined systems with M = 1 M⊙ (from left to right): {q = 1, Ω = 0°}, {q = 1, Ω = 90°}, {q = 0.5, Ω = 90°}, and {q = 0.2, Ω = 90°}. The top labels correspond to the period in terms of binary periods. The regions between the horizontal red lines correspond to real polar orbits (i.e. librating around i = 90°). Long-term integrations show that particles with initial conditions in the red/brown regions do not survive more than 105 periods. When Ω = 90°, stability islands appear at polar inclinations closer to the binary (between 2 and 3 aB). In the bottom panel, the colour scale shows the amplitude of libration of i. The continuous red lines are defined by Eq. (2) and correspond to the analytic limits between polar and prograde or retrograde orbits.

In the text

Fig. 6 Top panels: dynamical maps of inclined configurations with Ω = 90° with q = 1, 0.5, 0.2 (from left to right). The colour scale corresponds to the eccentricity variation Δe of the orbit. The continuous red lines are defined by Eq. (1) and correspond to the analytic limit between polar and prograde or retrograde orbits. Bottom panels: longitudinal cuts of the dynamical maps at i = 90°. The amplitude of Δe (Δi) is shown in red (blue). The grey shaded areas correspond to chaotic orbits with ⟨Y ⟩ ≫ 2.

In the text

	Fig. 7 Dynamical map for polar orbits for massless planets with initial a = 4.5 aB in the (eB, q) plane. The colour scale corresponds to Δe. The orbits identified with white dots are unstable. Chaotic orbits identified with the MEGNO indicator are indicated with grey squares.
In the text

Fig. 8 Dynamical maps of polar orbits (i = 90°, Ω = 90°) for massive planets in the (q, Mp) plane with eB = 0.5. The colour scale in the upper and lower panels corresponds to ΔΩ and Δ(Ω − ϖB), respectively. Initially, the planet is located at the equilibrium point i = 90°, Ω = 90°. The four initial conditions identified with white circles correspond to the orbits shown in Fig. 9. Each system was integrated for 20 000 yr (~ 650 000 binary periods).

In the text

Fig. 9 Evolution of the individual polar orbits of Fig. 8: (1) blue, (2) green, (3) brown, and (4) yellow. The secular coupling between eccentricity, inclination, and node is easily identified. The periods are of roughly 5 × 104 binary periods, depending on the couple of values (MP, q). The evolution of giant planets (MP > 10−4, i.e. brown region in Fig. 9) exhibits the circulation of node. According to the MEGNO indicator, orbits 1, 2, and 3 are regular (⟨Y ⟩~ 2); instead, 4 is chaotic after ~100 000 orbital periods of the binary (⟨Y ⟩ > 2).

In the text