Parametric study of polar configurations around binaries

TL;DR

This study investigates the stability and dynamics of polar circumbinary disks around binary stars using N-body simulations, revealing how resonance locations, planetary presence, and system specifics influence stability.

Contribution

It extends parametric analysis of polar configurations around binaries, identifying stability limits related to mass ratio, eccentricity, and resonance locations, with application to real systems.

Findings

Stability limits depend on mass ratio and eccentricity.

Giant planets affect the stability regions of polar orbits.

Polar configurations around HD 98800BaBb are influenced by nodal positions.

Abstract

Dynamical studies suggest that most of the circumbinary discs (CBDs) should be coplanar. However, under certain initial conditions, the CBD can evolve toward polar orientation. Here we extend the parametric study of polar configurations around detached close-in binaries through $N$-body simulations. For polar configurations around binaries with mass ratios $q$ below $0.7$, the nominal location of the mean motion resonance (MMR) $1~:~4$ predicts the limit of stability for $e_{B} > 0.1$. Alternatively, for $e_{B} < 0.1$ or $q \sim 1$, the nominal location of the MMR $1~:~3$ is the closest stable region. The presence of a} giant planet increases the region of forbidden polar configurations around low mass ratio binaries with eccentricities $e_B\sim0.4$ with respect to rocky earth-like planets. For equal mass stars, the eccentricity excitation $\Delta e$ of polar orbits smoothly increases with decreasing distance to the binary. For $q<1$, $\Delta e$ can reach values as high as $0.4$. Finally, we studied polar configurations around $HD~98800BaBb$ and show that the region of stability is strongly affected by the relative positions of the nodes. The most stable configurations in the system correspond to polar particles, which are not expected to survive on longer time-scales due to the presence of the external perturber HD~$98800AaAb$.

Figures

Tables

Table 1. Table 1: Dynamical configuration for the N 𝑁 N -body integrations of HD 98800.

\brSystem	$a$ (au)	$e$	$i$ (deg)	$\omega$ (deg)	$\Omega$ (deg)	Mass ($M_{\odot }$)
\mr$BaBb$	0.982	0.785	66.8	109.6	337.6	$Ba$=0.699, $Bb$= 0.582
Disc	2.6 - 4.6	0.0	26 or 154	0.0	17.0	- -
$AB$	54.58	0.52	88.6	64.	4.2	$Aa$+$Ab$=1.3
\br

Equations

$$F=\sqrt{\frac{1-e_{\mathrm{B}}^2}{1-5e_{\mathrm{B}}^2{\cos }^2\Omega +4e_{\mathrm{B}}^2}}.$$

$$\arcsin F<i<\pi -\arcsin F.$$

$$a_c=a_{\mathrm{B}}\left(1.6+5.1e_{\mathrm{B}}-2.22e_{\mathrm{B}}^2+4.12\mu -4.27e_{\mathrm{B}}\mu -5.09{\mu }^2+4.61{\mu }^2{e}_{\mathrm{B}}^2\right)$$

$$a_c$$

Full text

1 Introduction

Against all odds, the number of circumbinary planets is slowly increasing. To date, 22 circumbinary (P-type) planetary systems have been catalogued, with 11 transiting circumbinary planets detected by Kepler around nine binary star systems [1]. The binaries eccentricities range from quasi circular orbits (Kepler-47, $e_{\rm B}=0.02$ [2]) to highly eccentric (Kepler-34, $e_{\rm B}=0.52$ [3]). All the transiting binaries have short periods ($<30$ days). Moreover, the detected planets have almost circular coplanar orbits and semi-major axis ratios with the binary $a_{\rm p}/a_{\rm B}\sim 4$. Typical planets found around binary stellar systems have a radius of the order of $10$ Earth radii and orbital periods of about $160$ days (i.e. $a_{\rm p}\sim 0.35$ au) [4].

Circumbinary (P-type) planets are intrinsically difficult to observe through radial velocities, which translates in a discrimination in Doppler planet searches [5, 6]. Hence, these planets are more easily observed through transit methods — provided precise considerations on the geometry of transit are considered [7]. The occurrence rate of circumbinary planets is comparable to that of planets around single stars if the mutual inclination is always small ($<5\degree$), while this rate may be much larger if modest inclinations ($>5\degree$) are common [8].

For circumbinary planetesimals in coplanar orbits, eccentricities evolve on a dynamical timescale, which leads to orbital crossings even in the presence of gas drag. This makes the current locations of the circumbinary Kepler planets hostile to planetesimal accretion [9, 10]. Additionally, stellar-tidal evolution models of short-period binaries show that the binary orbital period increases with time [11, 12]. 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.

Planets form within circumbinary discs, that are typically assumed to be aligned with the binary [13, 14]. Misaligned circumbinary discs (CBDs) were considered as unlikely or at least transient; however, the observations of highly non-coplanar systems such as 99 Herculis [15], IRS 43 [16], GG Tau [17, 18], and HD 142527 [19] suggest otherwise. More importantly, the very first confirmation of a circumbinary gas-rich disc in a polar configuration HD 98800 [20] motivates our study. Previous works based on Smoothed Particle Hydrodynamics (SPH) simulations studied specific conditions for polar alignment of the disc around binaries [21, 22]. This mechanism has been further investigated providing an analytical framework to describe the polar alignment of CBDs [23, 24]. In particular, [25] conducted an extensive numerical exploration through SPH simulations. Lastly, the symmetry breaking between prograde and retrograde CBDs reported by [26, 27] could increase the likelihood of polar alignment.

For these reasons, it is important to further understand the stability of polar configurations around binary stars. In [27] we presented the most favourable conditions for regular movement of particles and planets around binaries. Here, we extend these results for a wider variety of configurations and analyse the case of HD 98800.

2 Dynamics around binaries

The secular evolution of a planet around a binary was presented in [28] and later investigated in [29], where they reported the coupling between the inclination and the node of particles. Then, [30] studied the circumbinary elliptical restricted and general three-body problem, giving an averaged quadrupolar Hamiltonian. In their formulation, the eccentricity of the outer particle/planet remains constant (i.e. $\Delta e=0$). This is not longer valid if the octupole expansion is considered or in full $N$-body integrations (as in [27] and this work).

An analysis of the stability of inclined massless particles was done by [31] where they consider different binary eccentricities and mass ratios between its components. 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 [32, 33, 34]. Interestingly, for near polar configurations, the inner binary can significantly excite the orbital eccentricity (up to about $0.3$ in some cases, as shown by [34]). Lastly, the evolution of polar orbits considering relativistic effects was studied by [35, 36]. The latter effects are beyond the scope of this work.

Here, we consider a binary system with total mass $M$ and individual masses $M_{1}$ and $M_{2}$, with a binary mass ratio $q=M_{2}/M_{1}$. To describe the motion of a P-type particle (planet) we use the Jacobi orbital elements semi-major axis $a$, eccentricity $e$, inclination $i$ (with respect to the binary orbital plane), mean longitude $\lambda$, longitude of pericentre $\varpi$ (alternatively argument of pericentre $\omega$), and longitude of the ascending node $\Omega$. The sub-index $B$ is used when referring to the binary orbit. Without loss of generality, we set initial conditions {$\lambda_{\rm B}=0\degree$, $\varpi_{\rm B}=0\degree$, $\Omega_{\rm B}=0\degree$}. The angles of a given particle are measured from the direction of the binary pericentre.

For a P-type particle the phase space has two equilibrium points at {$\Omega=\pm 90\degree,\,i=90\degree$}, where both angles librate. For any given pair of $(\Omega,\,i)$ it is possible to calculate the separatrix F as follows [30],

[TABLE]

Then, polar alignment is expected when

[TABLE]

Therefore, the higher the binary eccentricity $e_{\rm B}$ and the closer the node of the particle to $\Omega\sim 90\degree/270\degree$, the larger the region of polar orbits. In the next section, we numerically explore the stability for a wide range of $\{e_{\rm B},q,M_{\rm p}\}$ in polar configurations initially placed at $(i=90\degree,\Omega=90\degree)$.

3 Stability of polar orbits

Since polar configurations are a possible outcome for the evolution of misaligned gaseous discs around sub-au binaries [23, 24, 4, 27], it is relevant to identify which are the allowed configurations once the gas is dissipated, in a similar way as in [37] for the study of the coplanar case.

As mentioned previously, the binary sample of Kepler systems ranges from $0.2\leq q\leq 1$ and $0.05\leq e_{\rm B}\leq 0.52$. Therefore, we choose to explore the stability of massless polar particles for different binary mass ratios $q=\{0.2,0.5,0.7,1\}$.

In all our dynamical maps each orbit was integrated for at least 365 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 order to identify the dynamical regimes of movement, we calculated the Megno (Mean Exponential Growth of Nearby Orbits) value $\langle Y\rangle$ for each orbit [38], along with the maximum $\Delta\Omega$ and the maximum $\Delta e$ values attained during the dynamical evolution. We solve the $N$-body equations of motion numerically by means of a double-precision Bulirsch-Stoer integrator with tolerance $10^{-12}$.

In [27] we studied the behaviour of planets around a binary with $e_{B}=0.5$, and showed that planets with masses $M_{\rm p}>10^{-5}\,M_{\odot}$ exhibit rapid node circulation. In this case, the coupled oscillation with the eccentricity, leads the planet to describe a pulsating sphere around the binary. In Figure 1, we show the eccentricity excitation in the ($e_{B},q$) plane for different planetary masses ($1\,M_{\oplus}$, $17\,M_{\oplus}$, $1\,M_{\rm J}$), setting $a=4.5\,a_{\rm B}$. The value of $a_{\rm B}$ is similar to the observed configurations of Kepler systems $(a_{\rm B}\sim 0.1-0.2$ au). As expected, the chaotic orbits are correlated with high values of $\Delta e$. Moreover, we observe that the larger the planetary mass, the more significant the instability region. More specifically, Earth-like planets remain with low $\Delta\Omega$ and Neptune-like planets exhibit $\Delta\Omega\sim 90\degree$; whereas the node of Jupiter-Like planets always circulates.

In Figure 2, we show the eccentricity amplitude $\Delta e$ of polar orbits around a binary with $e_{\rm B}=0.5$, for different mass ratios $q$. This extends the region studied in [27]. The eccentricity excitation remains below $\Delta e<0.1$ only for binaries with $q\sim 1$. For different values of $q$, the particle experiences a significant eccentricity excitation, that strongly depends on its initial semi-major axis (contrary to the analytical model prediction in [30]). Higher values of $\Delta e$ eventually trigger instability, which leads to the disruption of the system. Hence, eccentricity excitation can have dramatic dynamical effects on long time-scales.

Finally, we choose to analyse the polar orbits in the plane ($a,e_{\rm B}$) for different mass ratios $q$. In Figure 3, we observe (almost) vertical stripes associated to the $1~{}:~{}N$ mean motion resonances (MMRs), where eccentricity is efficiently excited. For polar orbits, the $1~{}:~{}3$, $1~{}:~{}4$, and $1~{}:~{}5$ MMRs are the most prominent features. Their nominal location is displaced closer to the binary for increasing eccentricities $e_{\rm B}$. In Figure 4, we show a zoomed region of these maps for $q=0.2$ and $q=1$, where the limit of the stability region can be better appreciated. The latter is associated with the $1~{}:~{}N$ resonances.

The often-quoted study of [39] numerically derived a stability limit for circular coplanar particles around binaries:

[TABLE]

where $\mu=\frac{q}{1+q}$.

Alternatively, we found that in polar configurations with mass ratios below ($q<0.7$), the nominal location of the $1~{}:~{}4$ MMR predicts the limit of stability for $e_{\rm B}>0.1$; while, if $e_{\rm B}<0.1$ or $q\sim 1$, the nominal location of the $1~{}:~{}3$ MMR seems to be the innermost limit (see Fig. 3 and zoom in Fig. 4). When $e_{\rm B}\gtrsim 0.5$, the perturbations on the particles become more significant and only particles with $a\gtrsim 8\,a_{\rm B}$ remain unaffected by the perturbation. Then, using the MMR location, we can estimate a rule of thumb for the stability of polar orbits as:

[TABLE]

4 Inclined orbits around HD 98800 BaBb

The HD 98800 system is a hierarchical quadruple stellar system, composed by two pairs of binaries (called ‘$A$’ and ‘$B$’, or equally ‘$AaAb$’ and ‘$BaBb$’). The very eccentric binary $BaBb$ is well constrained [40] and hosts a bright circumbinary disc discovered by [41]. Using new data [20] derived the $AB$ orbit111We note that there are small systematic differences between the supplementary material and the main text.. The disc around ‘$BaBb$’ has an inner edge truncated by the binary $BaBb$, and the outer edge externally truncated by the companion $AaAb$. The orientation of the disc was initially thought to be coplanar with the inner binary, but higher-resolution observations suggest a different orientation recently confirmed by [20]. Two solutions were proposed for the disc  [20]: it is inclined either by $26\degree$ or by $154\degree$ with respect to the sky, which corresponds to an inclination with respect to the binary222The mutual inclination is obtained through spherical trigonometry: $i_{m}=\cos{i}\cos{i_{\rm B}}+\sin{i}\sin{i_{\rm B}}\cos{(\Omega-\Omega_{\rm B})}$ [42].. of $i_{m}=48\degree$ or $i_{m}=92.5\degree$ (respectively). A CBD in polar configuration is hence possible (and likely) in this system. The millimetric dust is observed to be on circular orbits with $2.5\text{ au}<a<4.6\text{ au}$ . Previous $N$-body simulations for the two possible inclinations showed that the dust is ejected in less than $1$ Myr, implying that that the dust observed with alma is embedded within a more massive gas disc that stabilises the disc against the exterior stellar perturbations of HD 98900$AaAb$ [20].

We study this system using initial conditions from Table 1, where the perturbing pair $AaAb$ is modelled as a single star like in [20]. We note that accurate orbit determination for binaries without a complete period is extremely difficult. As a matter of fact, orbital fits exhibit high correlation between the period and the eccentricity [43], depending on the phase covered by the observations. This correlation can be observed in the table given at http://www.ctio.noao.edu/~atokovin/stars/stars.php?cat=HD&number=98800.

In our previous section (see Fig. 3), the systems are initially located at the equilibrium configuration ($i=90^{o}$, $\Omega-\Omega_{B}=\pm 90^{o}$). The binary $HD~{}98800BaBb$ has a mass ratio $q\sim 0.83$, an eccentricity $e_{B}=0.785$, and $\Omega-\Omega_{\rm B}=-39.4^{o}$. Then, it is meaningful to analyse the stability of polar orbits for this kind of binary. In Figure 5, we show the stability of polar particles in the ($a$, $i$) plane without considering the external perturber $AaAb$, setting $\Omega-\Omega_{B}=\pm 90^{o}$ (left frame) and with the observational measured value of $\Omega-\Omega_{\rm B}=-39.4^{o}$ (right frame). In the latter case, the stable polar regions around the $1~{}:~{}4$ MMR disappear and only particles on polar orbits beyond the $1~{}:~{}6$ MMR are stable (that corresponds to $i=154\degree$). Our integrations show that even without the external perturber, the disc of dust is unstable for $a\lesssim 3.2$.

In Figure 6 we show the stability of polar particles around the binary HD 98800 BaBb in the ($a$, $e_{\rm B}$) plane, without considering the external perturber and setting $\Omega=337.6^{o}$ and $\Omega_{B}=17^{o}$ (i.e. $\Omega-\Omega_{B}=-39.4^{o}$). When the binary is eccentric ($e_{\rm B}=0.785$), the $1~{}:~{}6$ MMR location ($a\sim 3.5$ au) is the closest stable region. If the binary HD 98800 BaBb eventually circularises thought tidal interaction (see e.g. [44]), then the stability of closer polar orbits would be located around the $1~{}:~{}3$ MMR ($a\sim 2$ au).

Finally, in Fig. 7 we show the dynamical map in the ($a$, $i$) plane for the full four-body problem — coloured with the ejection time from the system — using the data from Table 1. The most stable conditions (i.e. with survival times of $\sim 1$ Myr) are found for almost polar configuration ($i=154\degree$), for $a\sim 4.2$ au and $a\sim 6$ au (beyond the measured location of the disc of dust). Interestingly, if new orbital determinations put the exterior perturber on a wider orbit, then the polar orbits around the binary $BaBb$ would become stable on longer timescales.

5 Conclusions

Here we extended the previous study by [27] on polar configurations around eccentric binaries. The stability of single polar planets around binaries is weakly dependent on the planetary mass (see Fig. 1). There is however a planetary mass threshold ($M_{\rm p}\gtrsim 10^{-5}M_{\odot}$) above which the planet exhibits rapid node circulation (as shown in [27]).

We find that the eccentricity excitation strongly depends on the distance to the binary. Then, assuming planets form in the external disc regions and then migrate inwards, this excitation eventually triggers instability. Alternatively, the eccentricity excitation remains very low for $e_{\rm B}<0.4$ when $q<0.7$. Generally, for massless particles, the orbits with $0.2<e_{\rm B}<0.4$ exhibit chaotic Megno (see Fig. 2).

For polar orbits around binaries with mass ratios $q$ below $0.7$, the nominal location of the $1:4$ mean motion resonance (MMR) predicts the limit of stability for $e_{\rm B}>0.1$; while if $e_{\rm B}<0.1$ or $q\sim 1$, the nominal location of the $1~{}:~{}3$ MMR is the closest stable region to the binary (see Fig. 3).

Finally, we analysed the HD 98800 system where there is a circumbinary disc around $BaBb$ being currently perturbed by an outer binary $AaAb$. In particular, we studied the stability of the system with and without the additional binary, for any arbitrary inclination. We find that the near polar configuration for the circumbinary disc around $BaBb$ is the more stable among the all possible disc inclinations. We also conclude that — in the absence of gas — the instability would not be due to the polar configuration itself, but rather to the gravitational perturbations of the outer binary $AaAb$ (see Fig. 7).

\ack

This work has been supported by research grants from CONICET and Secyt-UNC. $N$-body computations were performed at Mulatona Cluster from CCAD-UNC, which is part of SNCAD-MinCyT, Argentina. NC acknowledges financial support provided by FONDECYT grant 3170680 and from CONICYT project Basal AFB-170002. 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.