© The Authors 2023

1 Introduction

Circumbinary planets (CBPs) discovered around transiting binaries challenge our knowledge of planet formation and migration in the disc around binaries. There are 13 transiting CBPs orbiting 11 Kepler-eclipsing binaries. Recently, the first circumbinary planet was reported using data from the TESS mission (TOI 1338, Kostov et al. 2020). Most of the planets lie near their instability boundary at about three to five binary separations in regions surrounded by strong resonances. Recent studies argue that the occurrence rate of giant Kepler-like CBPs is comparable to that of giant planets in single-star systems (see Martin et al. 2019, and references therein).

Transiting binaries that host planets have short periods (<30 days) (Schwarz et al. 2016) with eccentricities ranging from quasicircular (such as Kepler-47, Orosz et al. 2012b) to highly eccentric orbits (such as Kepler-34, Welsh et al. 2012). Typical planets detected around binary stellar systems have a radius of about 10 Earth radii and orbital periods of about 160 days on almost circular coplanar orbits (i.e. a_p ~ 0.35 AU). Radius, period, and coplanarity values are affected by observational biases (Martin 2018). However, long-period CBPs could very well exist, with various masses and sizes.

There is some evidence that the N/1 mean motion resonances (MMRs) are related to the planetary parking location and their evolution in the disc, as they lead to strong regions of chaos around the binaries (Gallardo et al. 2021). For example, Zoppetti et al. (2018) applied a simplified model to Kepler-38 and suggested a capture in 5/1 MMR with subsequent tidal evolution of the planet outside MMRs. Other authors have performed fine-tuning hydrodynamical simulations, placing some constraints on the protoplanetary discs that allow migration of CBPs, but these simulations do not show capture at MMRs (see e.g. Thun & Kley 2018; Penzlin et al. 2019, 2021, and references therein). In all these simulations, the authors find that the final parking locations of CBPs depend on the binary orbital parameters, which affect the CBD morphology. Essentially, the planets migrate to the inner edge of the disc cavity and remain there. But this is a problem because the parking locations seem to be too large when compared to recent observations. In their pioneering work, Nelson (2003) showed a temporary capture in MMR 4/1 around a binary. However, the planet was placed inside the disc gap and the system was integrated for a few hundred binary orbital periods only.

In situ planet formation at these close orbits around binaries is extremely unlikely. Strong perturbation forces close to the binary inhibit planetesimal and dust accretion (Moriwaki & Nakagawa 2004; Paardekooper et al. 2012; Silsbee & Rafikov 2015; Meschiari 2012), and turbulence induced by hydrodynamical parametric instabilities can dramatically reduce pebble accretion efficiencies (Pierens et al. 2020). A more likely scenario is planet formation in the outer disc and subsequent disc-driven migration to the observed close orbits (Pierens & Nelson 2008; Bromley & Kenyon 2015). However, this assumption does not explain the particular stopping positions of most CBPS.

In planet-forming discs, gas is partially supported against gravity by its pressure gradient. Therefore, the surrounding planetesimals experience aerodynamic drag, which translates into eccentricity damping. This mechanism is thought to increase their stability around single stars (Kley & Nelson 2012; Sutherland & Kratter 2019). However, the disc surrounding a binary may have a non-axisymmetrical distribution and become eccentric (Kley et al. 2008; Ragusa et al. 2017; Thun et al. 2017; Poblete et al. 2019; Hirsh et al. 2020; Paardekooper et al. 2022) because of the effects of the secondary star.

Several studies have already investigated planet migration in circumbinary discs (Nelson 2003; Pierens & Nelson 2008, 2013; Kley & Haghighipour 2014, 2015; Kley et al. 2019). Because of gravitational torques, the binary clears out the inner cavity and planet migration is stopped at the inner edge of the disc at a location close to the binary but not yet matching the observations in most cases studied. Numerical simulations show that the inner disc shape and therefore the final planet orbital parameters depend on disc parameters such as viscosity and scale height, and the mass ratio and eccentricity of the binary (Mutter et al. 2017; Thun & Kley 2018; Ragusa et al. 2020)

Our aim is to study planetary migration around different kinds of binaries and to elucidate the possible physical parameters of discs that can explain the final parking position of the planets. Taking into account the two exoplanets recently presented by Esmer et al. (2022), 42 circumbinary exoplanets have so far been discovered¹. We mainly rely on the findings of the Kepler transit survey when examining trends of circumbinary planets. This is because, in contrast to the many limitations of the suggested eclipse time variation planets, it is the only sample that is sufficiently large for preliminary population studies and also contains reliable discoveries. Furthermore, by restricting ourselves to a single observation method, only a single observation bias needs to be taken into consideration (Martin 2018). In the remainder of this paper, we restrict our study to CBPs detected by transit methods only.

Table 1 shows the planetary and star parameters from known transiting CBPs. We observe that typical planets around binaries are quite circular (e_mean ≃ 0.06), which could be a consequence of the interaction with their natal protoplanetary disc. All the planets have radii of greater than 3 R_⊕. Planet masses are typically poorly constrained or compatible with zero. We take metallicities from discovery papers, because we find that the Exoplanet.eu database² has missing values. For instance, the discovery paper of TOI-1338 (Kostov et al. 2020) reports [Fe/H] = −0.01 ± 0.05 dex (measured with HARPS) while its metallicity in the Exoplanet.eu database is [Fe/H] = −0.4 ± 0.1 dex.

Figure 1 shows the masses of exoplanets presented in the catalogue of Otegi et al. (2020, hereafter OC2020) as a function of their radius (grey markers). This catalogue contains data for planets with masses lower than 120 M_⊕, with reliable and robust mass and radius measurements. Known CBPs are also included (red markers) for comparison. The distribution of radii of exo-planets is bimodal (David et al. 2021), with peaks at ~2 R_⊕ and ~13 R_⊕, while the CBP distribution has a peak close to the valley of this distribution, at ~7 R_⊕. On the other hand, although the mass distribution of exoplanets can also be categorised as bimodal with peaks at ~9 M_⊕ and ~260 M_⊕, the mass distribution of CBPs is harder to categorise because of the large errors in their estimations. However, this distribution exhibits a peak close to the valley of the exoplanet mass distribution, at ~60 M_⊕. From these comparisons it is possible to reaffirm that the outcome of the formation and evolution of CBPs could be different from that of exoplanets around single stars. It is worth mentioning that Kepler-47 is the only multi-planetary circumbinary system, and Kepler-64 (also known as PH1 b) is a circumbinary exoplanet in a four-star system. The right panel in Fig. 1 shows a zoom onto the region around the CBPs. All CBPs reside above the theoretical H₂O composition line, and therefore could be considered as volatile-rich. These planets may therefore have formed beyond the disc ice line at large semi-major axis before migrating towards the inner disc regions. From this figure, it would appear that binaries with larger eccentricities tend to host larger and more massive planets. A similar correlation seems to be fulfilled when one considers the binary mass instead. By way of comparison, this figure also shows a fit for the mass–radius relationship (valid for volatile-rich planets) presented in OC2020, given by:

$$\left(\frac{R}{R_{\oplus }}\right)=\left(0.7\pm 0.11\right){\left(\frac{m}{M_{\oplus }}\right)}^{\left(0.63\pm 0.04\right)},$$(1)

where R is the planet radius and m its mass. Although it is evident that the CBPs do not follow the trend of this fit (red curve in right panel in Fig. 1), it is difficult to calculate a reliable fitting function because of the paucity of available data.

Figure 2 shows the radii of known CBPs as a function of their period ratio (P/P_B ≡ n_B/n) with their host binaries. None of the observed CBPs have a period ratio of below 5, and most of them have a period ratio of between 5 and 7. It also seems to be the case that the larger the binary eccentricity, the larger the period ratio of the planet. These characteristics seem to be independent of the planet mass. In the following section, we describe our numerical simulations to understand planetary migration around binaries.

Fig. 1 Mass–radius relation of transiting circumbinary exoplanets and circumstellar exoplanets. Left panel: radii of exoplanets presented in OC2020 (grey) and CBPs (red/black) as a function of their mass. Their respective error estimations are presented with error bars. In the case of K38, K47b, K413, and K453, we mark their maximum mass reference value with cross symbols (×). The light-blue, orange, and green lines denote the theoretical H2O, Fe, and Earth-like composition lines. The curves outside the box represent the kernel density estimation of the mass (top) and radius (right) of the scattered data points. The data bandwidth for density estimation is calculated with Scott’s Rule (Scott 2015) using the scipy Python package. The vertical (top) and horizontal (right) dashed lines denote the locations of the maximum values in each distribution. All distributions are normalised to their own total counts. Right panel: zoom onto the CBP region shown in the left panel. The exoplanets presented in OC2020 are shown as grey circles, while the CBPs are coloured according to their binary mass. The size of each CBP circle is proportional to the eccentricity of the binary. The red shaded area denotes the mass–radius relation (the red solid line being its reference values) presented in OC2020

Fig. 2 Period ratio of known circumbinary planets as a function of their radii. The size of each point represents the planet mass, and the colour scale indicates its binary eccentricity. Horizontal grey lines denote the location of the resonances N/1. The plot is truncated to P/PB|imax = 11.

2 Methods

Based on N-body integrations without dissipation forces, we construct dynamical maps in the plane a–e to highlight the location of MMRs and the process of capture when performing the long-term N-body with Stokes-like dissipation forces. We consider a binary system with masses m₀ and m₁, respectively, and semi-major axis a_B = 0.1 AU. The mass of the planet orbiting the binary is equal to m = 93 M_⊕. We choose a Jacobi system for the coordinates. The orbital elements a, e, M, and ω are the semi-major axis, eccentricity, mean anomaly, and argument of the pericentre, respectively. We use the subscript B to refer to the binary orbit. The binary and planet orbits are coplanar.

The N-body simulations were performed with a Burlisch–Stoer integrator with adaptive step size, which was modified to independently monitor the error in each variable; this imposes a relative precision of better than 10⁻¹³. Integrations are stopped when the distance from the planet to any star is less than the sum of their radii, or when the planet is ejected from the system after scattering between bodies (>2 AU). To construct the dynamical maps (without considering external forces), the total integration time was set to encompass several periods of the orbital secular variations (i.e. 600 yr ≈ 20 000 binary periods).

The initial conditions for the circumbinary dynamical maps are coplanar, with mean anomaly M_B = M = 0° and argument of pericentre ω_B = ω = 0°. To better sample the structure of the resonances, we calculate ∆e as the amplitude of maximum variation of the orbital eccentricity of the planets during integration, with ∆e = (e_max − e_min). This indicator correlates with the chaos indicator MEGNO (Cincotta & Simó 2000), and we verified that higher values of ∆e (>0.4) correspond to chaotic orbits (i.e. MEGNO >2).

In the dynamical maps, the locations of N/1 resonances were calculated using the method developed by Gallardo et al. (2021), using a semi-analytical model for planetary MMRs around binary stars. The code used can be found online³. We assume that the planet is aligned with the binary and that the initial mean anomalies of the binary and the planet are equal to zero. We recall that, as shown by Gallardo et al. (2021), the width of the resonance for coplanar orbits depends on the initial relative value of δϖ = ϖ_B − ϖ = ω_B − ω. Given that the pericentre of the planet precesses with a short period, the dynamical maps are similar for any initial value of δϖ as we approach the binary. It is also worth mentioning that for resonances closer to the binary, the nominal location of a given MMR depends on the mean anomaly (see e.g. Fig. A.2. in Giuppone et al. 2022).

To simulate the disc-induced migration of the planets, we added an ad hoc external force to an N-body integrator. This additional term behaves as a Stokes non-conservative force:

$$\frac{{\text{d}}^{2}r}{\text{d}{t}^2}=-C\left(\upsilon -\alpha {\upsilon }_{\text{c}}\right),$$(2)

where r is the position vector of the planet (in a Jacobian reference plane), υ its velocity vector, and υ_c the circular velocity vector at the same point (Beaugé et al. 2006). With this formulation, the coefficients C and α are defined as

$$\begin{array}{cc}C=\frac{1}{2{\tau }_a}+\frac{1}{{\tau }_e},& \alpha =\frac{2{\tau }_a}{2{\tau }_a+{\tau }_e},\end{array}$$(3)

where τ_a and τ_e are the characteristic semi-major axis and eccentricity damping timescales, respectively.

Following Beauge & Ferraz-Mello (1993), at first order in eccentricity and for a single planet, the effects of the previous force in the semi-major axis and eccentricity of the body can be described as:

$$\begin{array}{cc}a\left(t\right)=a_0\exp \left(-\frac{t}{{\tau }_a}\right),& e\left(t\right)=e_0\exp \left(-\frac{t}{{\tau }_e}\right),\end{array}$$(4)

where a₀ and e₀ are the conditions at the beginning of the integration. The a-folding and e-folding times for a and e are denoted |τ_a| and |τ_e|, respectively, and have the following expressions:

$$\begin{array}{cc}{\tau }_a^{-1}=2C\left(1-\alpha \right),& {\tau }_e^{-1}=C\alpha \end{array}.$$(5)

Following this formalism, the accelerations from tides were incorporated into our N-body code, as in Ronco et al. (2020) for other external forces.

Some authors set the initial parameters τ_a, τ_e from the binary and disc properties (Zoppetti et al. 2018; Secunda et al. 2019, 2020; Martin & Fitzmaurice 2022), and then evolve them through simulations according to some defined relation (Cresswell & Nelson 2008). These relations were obtained mainly from hydrodynamical simulations of discs around single stars, and their application to circumbinary planets should be done with caution (Chrenko et al. 2018). In order to parameterise τ_e through simulations, some authors use the simple relation

$${\tau }_e=\frac{{\tau }_a}{K},$$(6)

where K is a constant with a typical value K = 10 (Lee & Peale 2002; Rein 2012; Martin & Fitzmaurice 2022).

In our work, we extend the parametric study presented by Martin & Fitzmaurice (2022), increasing the range of explored parameters. This approach is mainly motivated by the fact that the torque equations around binary systems have not yet been fully resolved; especially near the binary baricentre where important MMRs like 5/1, 4/1, and 3/1 also play an important role in the dynamics. Instead, for our simulations, we set the initial values of τ_a and τ_e independently of each other, and they remain constant throughout all simulations.

We vary the following parameters in our simulations: the binary eccentricity and mass ratio (e_B, q), the mass and initial eccentricity of the planet (m, e), and the characteristic timescales of the semi-major axis and eccentricity dampings (τ_a, τ_e). In all our simulations, the initial semi-major axis of the planet is set at a = 1.5 AU, which is sufficiently far from the binary star, which interacts both with the planet and the disc.

Table 2 summarises the different values adopted for each parameter in our simulations. Due to the fact that the expected lifetime of each simulation is roughly min (τ_a, τ_e), we set the output times for each simulation as follows:

$$\text{d}{t}_{\text{out}}=0.5\times {10}^{{\log }_{10}\left(\min \left({\tau }_e,{\tau }_a\right)/\text{yr}\right)-2}\text{yr.}$$(7)

With this configuration, we are able to avoid storing large amounts of output data without subsampling.

Instead of shutting down the external ad hoc force (which would take into account the dissipation of the disc), we kept it active throughout the run. By doing so, the disc does not dissipate. We choose this configuration because uncertainty remains over the estimated lifetime of circumbinary discs (Alexander 2012; Shadmehri et al. 2018; Ronco et al. 2021), although disc lifetime measured observationally is about 3 × 10⁶ yr (Ribas et al. 2014).

3 Dynamical maps and drag of the planets

Figure 3 shows six dynamical maps in the plane (a, e) around a binary with mass ratio q = {0.2, 1} (top and bottom frames, respectively) and eccentricities e_B = {0.05, 0.3, 0.5} (left, centre, and right frames, respectively). Each dynamical map consists of a grid of 300 × 100 initial conditions integrated for 600 yr. The colour scale varies with ∆e, increasing from blue (regular motion) to red (chaotic motion). This figure also shows the location of the most important N/1 MMRs, and the stability limit

$$\begin{array}{lll}{a}_c\hfill & =\hfill & a_{\text{B}}({1.60}_{-0.04}^{+0.04}+{5.10}_{-0.05}^{+0.05}{e}_{\text{B}}-{2.22}_{-0.11}^{+0.11}{e}_{\text{B}}^2\hfill \\ \hfill & \hfill & +{4.12}_{-0.09}^{+0.09}\mu -{4.277}_{-0.17}^{+0.17}\mu e_{\text{B}}-{5.09}_{-0.11}^{+0.11}{\mu }^2\hfill \\ \hfill & \hfill & +{4.61}_{-0.36}^{+0.36}{\mu }^2{e}_{\text{B}}^2),\hfill \end{array}$$(8)

from Holman & Wiegert (1999), where μ = m₁/(m₀ + m₁), and using the critical eccentricity approximation,

$$e_{\text{c}}=0.8\left(1-\frac{a_c}{a}\right),$$(9)

from Quarles et al. (2018). These widely used empirical criteria (Yamanaka & Sasaki 2019; Zoppetti et al. 2020; Gallardo et al. 2021; Martin & Fitzmaurice 2022) roughly define the smallest stable orbit around a binary system.

In all six maps, it is possible to identify a minimum eccentricity variation region (secular mode, Moriwaki & Nakagawa 2004; Paardekooper et al. 2012; Zoppetti et al. 2019). The location of the ‘capture’ mean eccentricity of this region can be approximated by

$${\langle e\rangle }_{\text{cap}}=\sqrt{e_{\text{f}}^2+{\in }^2\langle \frac{T}{L^{*}}\rangle },$$(10)

where

$$e_{\text{f}}=\frac{5a_{\text{B}}{e}_{\text{B}}\left(3e_{\text{B}}^2+4\right)\left(m_0-m_1\right)}{8a\left(3e_{\text{B}}^2+2\right)\left(m_0+m_1\right)}$$(11)

is the forced eccentricity, and ϵ²〈T/L^*〉 is an eccentricity term that can be approximated with the following analytical expression (Paardekooper et al. 2012; Zoppetti et al. 2019):

$${\in }^2\langle T/L^{*}\rangle =\frac{9}{16}\frac{m_0^2{m}_1^2}{{\left(m_0+m_1\right)}^4}{\left(\frac{a_{\text{B}}}{a}\right)}^4\left(1+\frac{34}{3}{e}_{\text{B}}^2\right).$$(12)

The values for Eq. (10) are shown in solid magenta lines in Fig. 3.

Even though the Holman & Wiegert (1999) limit considers the growth of the chaotic region (white area of the maps) for binaries with higher eccentricity or mass ratio, it fails to properly enclose the stable region. For example, in the upper right frame of Fig. 3, there is a roughly stable region (of the secular branch) when a ϵ [0.36, 0.42] AU and e ~ 0.1, which is considered unstable by this limit. From Fig. 3, it is also important to note the presence of isolated regions that are not fully unstable (∆e ≲ 0.3) near a ϵ [0.35, 0.4] AU and a ϵ [0.33, 0.35] AU in the upper and lower right frames (respectively), also incorrectly catalogued as unstable by this limit.

This figure also shows the migration trajectory of a planet of mass m = 5 M_⊕ with initial eccentricity e = 0.01, τ_a = 10⁵ yr, and τ_e = 10⁶ yr. The trajectories are plotted back to the last non-ejected state, and in all six simulations the planet got ejected at ~1.5 × 10⁵yr. In addition to showing every (a_i, e_i) point for a planet, its simple moving average (SMA) is also shown as a red line. This estimator is defined as:

$$SMA{\left(X,N\right)}_j=\frac{{\displaystyle {\sum }_{i=j}^{j+N}{X}_i}}{N},$$(13)

where N is the SMA window. N = 20 for every SMA calculation in this work.

We generally observe that the planet migrates inwards through the secular mode, until some MMR captures it. When capture occurs, the planet eccentricity increases with a small variation in the semi-major axis. For the higher binary eccentricity, or lower binary mass ratio, the secular mode is located further away from the binary and has higher eccentricity. This relationship causes the capture of exoplanets at higher eccentricities by higher order MMRs in low q or high e_B binaries. It is worth noting that the planets in the upper and lower right frames of Fig. 3 were able to reach the stable islands mentioned previously, migrating through the unstable region in between. Eventually, these planets are ejected (see Appendix A).

To illustrate the dependence on mass and initial eccentricity of the planet, we show two examples for the same binary (q = 0.2, e_B = 0.05) and dissipation parameters (τ_a = 10⁶ yr, τ_e = 10⁷ yr), but setting different initial conditions for the planet. The upper frame of Fig. 4 shows the final evolution of three planets with mass m = 5 M_⊕ and initial eccentricities e = 0.001, 0.01, and 0.1. In this case, the planet with eccentricity e = 0.1 is not able to damp its eccentricity down to the secular branch, and therefore it interacts with the MMR 6/1 within a region of non-negligible resonance width. As a consequence, the planet is captured by this resonance and subsequently ejected.

On the other hand, planets with eccentricities e = 0.001 and e = 0.01 migrate with a low average eccentricity. Even though both of them migrate through the secular branch, the first is captured in MMR 4/1, while the other in MMR 5/1. These results suggest that there is a relationship between the initial eccentricity of the planet and the resonance that captures it. More specifically, it appears that exoplanets with higher initial eccentricity are captured by higher order resonances. We verify that the MMR angle librates when capture occurs. This is shown in Appendix B for these three specific examples.

The lower frame of Fig. 4 shows the final evolution of three planets with eccentricity e = 0.01 and masses m = 5 M_⊕, 20 M_⊕, and 93 M_⊕. In this case, the mass of the planet also affects the resonance capture. Although it might appear that the results shown in this panel are suggesting that lower mass planets are ejected by higher order resonances, a simple correlation between capture and mass is not evident.

It is not practical to show all our results in this kind of figure because, as mentioned, we explored 4200 migration simulations. Therefore, in the following section we take a statistical approach to identifying the escape resonance.

Fig. 3 Dynamical maps in the plane (a, e), where each initial condition is integrated for 600 yr and the colour scale corresponds to the ∆e indicator. Solid black lines denote the width of the most prominent N/1 MMR and the dashed vertical lines indicate their nominal location, both calculated using Gallardo et al. (2021). For each of these maps, the binary has q = 0.2 (top frames) and q = 1 (bottom frames), and eB = 0.05 (left frames), eB = 0.3 (centre frames), and eB = 0.5 (right frames). The grey shaded area denotes the stability limit (the blue dot-dashed line being the reference value) calculated by Holman & Wiegert (1999), with its critical eccentricity approximated by Quarles et al. (2018). The solid magenta lines denote the approximation of the ‘capture’ mean eccentricity (Zoppetti et al. 2019). The light-blue dotted line denotes the simulated (a, e) evolution of a planet with m = 5 M⊕ and e = 0.01, setting τa = 105 yr and τe = 106 yr. The red solid line denotes the simple moving average of the planet trajectory, with a window of 20 points (dtw = 105 yr).

Fig. 4 Moving average trajectory in the plane (a, e) of three individual simulations with the same τa = 105 yr, τe = 106 yr, q = 0.2, m = 5 M⊕, and eB = 0.05, but different initial eccentricity. Vertical dashed lines denote the location of the 4/1, 5/1, and 6/1 MMRs. The solid magenta lines denote the approximation of the ‘capture mean eccentricity (Zoppetti et al. 2019).

Fig. 5 Distribution of R values created with KDE, and a histogram of their closest MMR value (Rint). The values above each Rint bin indicate the fraction of the total 4200 simulations found in that bin as a percentage. The R value KDE distribution is re-scaled for improved visualisation.

4 Final fate of migrating circumbinary planets

4.1 Kernel density estimation for ejecting MMRs

As shown above, the value of the MMR that captures the migrating planet in each simulation is the most relevant parameter for predicting/understanding its final fate. Given that the ad hoc Stokes force introduced in our simulations is never turned off, it is not surprising to find that the vast majority of exoplanets are eventually ejected by some MMR. Using a more realistic migration prescription where the Stokes-like force eventually vanishes, we expect planets to remain trapped in those MMRs without being ejected.

Because of the large number of simulations, it is impossible to inspect each simulation carefully to retrieve the MMR that captures and ejects the planet. Therefore, we designed a simple method to calculate this from the trajectory data:

1. Calculate the instantaneous semi-major axis ratio R = (a/a_B)^1.5 for each output time, which is equivalent to the period ratio P/P_B.

2. Define a Gaussian kernel density estimation (KDE) from the calculated R distribution, with a fixed bandwidth (=0.03, in this work).

3. Calculate the KDE value for fixed bins (points) and retrieve the bin with the highest score.

The R bins (note that this variable is dimensionless) used for these calculations consist of 250 points between 2.5 and 10.5 (both included), with a step of ~0.032. This binning contains all possible values of R for our simulations, and allows calculation of the proximity between each R value and its closest MMR. To avoid taking into account the large number of high values of R (corresponding to the first interval of planet migration), we set an upper bound on the semi-major axis for the calculation of the value R at a_max = 1 AU.

Multiple results of this method were compared with visual estimates of the trajectories, showing very close agreement. However, because this method is restricted by the dt_out of each simulation, it can generate an incorrect classification of R in simulations that do not retain their planet captured in MMR for timescales longer than dt_out.

Figure 5 shows the distribution of the values R obtained using this method (renormalised for proper visualisation) and the distribution of their closest MMR value. The numbers at the top of each bin indicate the percentage of ejected planets for a given R_int. Almost all values of R are slightly shifted to the right of their respective MMR nominal location. This result arises from the fact that almost all planets first encounter their ejecting MMR from the right (because they always migrate inward). From this figure, we can confirm that this distribution is far from uniform, and the MMR 4/1 is the one with most planet captures and ejections. Of the 4200 planets simulated, only 346 (8.2%) reached the MMR 3/1, where they were ejected. On the other hand, 1553 of the planets (37%) were ejected at the MMR 4/1. For R ≥ 5, the number of planets ejected decreases with lower MMR values, starting with the minimum value of MMR 9/1, which manages to eject only 12 planets (<1%).

It is worth mentioning that some of the planets with final R values close to the centre between two adjacent MMRs were probably ejected by the contribution of these two MMRs. For example, several R ~ 3.5 values correspond to planets ejected by the contribution of MMRs 3/1 and 4/1. Generally, we observe that the process of resonant capture is responsible for planetary ejection.

To establish the final fate of the planets in all simulations in Table 2, the relevance of each varied parameter must be understood and categorised. In the following section, we analyse the relationship between the parameters and the resulting R value for each simulation.

4.2 Relevant parameters for predicting capture by MMRs

One method to evaluate the importance of each feature⁴ within the simulation dataset is to analyse the index of dispersion (D) of a target over all single values of each feature (Cox & Lewis 1966). This index measures how closely a set of observed occurrences resembles a typical statistical model in terms of clustering or dispersion. This index is defined as the ratio of the variance σ² over the mean value μ: D = σ²/μ. In our case, we take the instantaneous semi-major axis ratio (R) as the target value.

Table 3 shows the minimum, mean, and maximum index of dispersion of the R values, calculated among the indices obtained for each single value of each parameter in Table 2. τ_a/τ_e is also considered as an additional (derived) parameter. The distribution of R values for each e_B value has the lowest dispersion index of all the features. Therefore, e_B is the parameter that is the most closely associated with the final value R of each planet. This information is connected to the fact that the binary eccentricity is directly related to the position of the innermost stable MMR, as mentioned in Sect. 3 (Eqs. (8)–(12)). The remaining features have indices that are comparable to each other, and so it is not possible to clearly determine the order of their respective importance with this method.

To visualise the effect of varying each parameter of Table 2, we present average resonance tables in Fig. 6. This figure shows the effect of different initial e (top frame), q (second frame), m (third frame), and log₁₀(τ_a/τ_e) (bottom frame) on the average R value (${\mu }_{\text{R}}=\overline{R}$) and its dispersion (σ_R). All grids have e_B as their primary axis because of its strong relationship with R.

From the top panel, we can determine that planets with higher e (= 0.1) tend to be captured by slightly higher order MMRs (ranging from 4.5 to 7.4). This fact is related to the resonant capture of planets with high eccentricity, which is discussed at the end of Sect. 3. There are no appreciable differences in the values of μ_R between planets with low eccentricity (e < 0.1).

From the second panel, it is possible to infer that the dependence of μ_R on q is weakly related to the binary eccentricity. For a large value of the mass ratio, the increase rate of μ_R with respect to e_B has a smaller slope compared to systems with a smaller mass ratio. In this case, the slope for simulations with q = 1 is ~(5.27 ± 0.48)μ_R/e_B, while that for simulations with q = 0.2 is ~(6.49 ± 0.60)μ_R/e_B.

From the third panel, the value of μ_R does not vary as the mass of the planet changes. Therefore, regarding the capture in MMRs, the planetary mass is a parameter of minor significance. From the lower panel it is possible to determine the relationship between log₁₀(τ_a/τ_e) and μ_R. The value of μ_R decreases as the value of |log₁₀(τ_a/τ_e)| increases, regardless of whether τ_a is higher or lower than τ_e. This effect is more noticeable as the binary eccentricity increases. It is interesting to note that, although the values of μ_R are similar, their dispersion values σ_R are slightly higher for the simulation groups with τ_a ≤ τ_e.

4.3 Methods to predict the MMR at which ejection occurs

The large number of simulations performed allowed us to carry out a statistical analysis to determine whether it is feasible to predict the MMR at which a CBP will be caught, depending on the initial values of the parameters examined. For this purpose, we implemented a simple random forest (RF) model.

As an ensemble learner, a RF generates a large number of classifiers and aggregates their results. This algorithm builds several classification (or regression) trees, each of which will be trained on a bootstrap sample of the initial training data and will search on a randomly chosen subset of input variables to evaluate the splitting (Breiman 2001). For a classification task, the class chosen by the largest number of trees is the RF output. To predict a variable, this predictive modelling tool pulls information from input data sets and seeks to identify new associations (Ivezic et al. 2014). In this work, the closest nominal MMR of an exoplanet is predicted by the classification task of this tool given initial condition parameters as input. For this task, we used the algorithm RandomForestClassifier⁵ from the scikit-learn Python package.

From the sample of 4200 simulations, we designed a training set made up of 75% of them (3150) and a test set made up of the remaining 25% (1050). The data columns (input parameters) of each set are [e_B, e, τ_a, τ_e, q, m, τ_a/τ_e], and their target is the integer closest to each value R, named R_int ≡ MMR. For this simple random forest classifier (RFC) model, we used the default hyper parameters from scikit-learn.

Figure 7 and Table 4 show the confusion matrix and a classification report (respectively) of our classifier applied to the test data set. This data set was never seen by the respective estimator during training. The metrics presented here are precision (TP/(TP + FP)), recall (TP/(TP + FN)), and F1 score (2 × precision × recall/(precision + recall)), where T(F)P(N) are the number of true (false) positive (negative) values. These metrics are some of the most widely used to measure the general performance of a trained model (Bengfort & Bilbro 2019; Schlecker et al. 2021; Audenaert et al. 2021). Table 4 also shows the support for each class, which is the actual number of instances of each class in the test data set.

The precision rate, which can be seen as a measure of the exactness of a classifier, is greater than 70% for all classes (except for R_int = 7, which has a rate equal to 69%). Moreover, the accuracy score of the entire set, calculated as the weighted average (averaging the support-weighted mean per label) of the precision rates, is 79%. Although this metric indicates how often the model makes a correct prediction, given the severely unbalanced nature of the data set, it is insufficient as a performance indicator. Therefore, we assess our model’s performance and examine the resulting different kinds of errors.

We can see in Fig. 7 that there is a significant difference between the colors of the diagonal and those of the off-diagonal entries, which indicates that our model has good performance. To retrieve the completeness of the classifier, which can be described as its ability to correctly classify all true instances of a given class, the recall rate is measured. Table 4 shows that, except for R_int = 6 and R_int = 9, all classes have a recall rate of greater than 70%, with 89% being the highest. The highest misclassification rate (33%) occurred between classes 8 and 9 (as one of the only three instances of R_int = 9 was classified as R_int = 8), and the highest accuracy rate (91%) was obtained for class 3. It is interesting to note that of all misclassifications, the vast majority occurred between adjacent classes (e.g. R_int = 4 with R_int = 5, or R_int = 8 with R_int = 7). This result indicates that the classifier was able to generate close relationships between the parameter set and the target classes.

Using this model, it is not possible to establish a parametric relationship between R_int and the input parameters (Ulmer-Moll et al. 2019). However, we are able to measure the importance of each feature of the data set according to the classifier. ‘Permutation feature importance’ is a technique for inspecting any fitted estimator when the data are tabular. This method retrieves the drop in model performance caused by randomly shuffling a single feature value (Breiman 2001). This procedure demonstrates the extent to which the model depends on the feature by breaking the feature–target relationship (Lu et al. 2020). The benefit of this method is that it can be calculated numerous times with various feature permutations and is model-independent. Figure 8 shows the importance of trained RFC characteristics, which can be understood as the contribution of each individual parameter to the final prediction. These values were calculated using permutation feature importance, with 15 permutations.

With almost 52% importance, the eccentricity of the binary is the most sensitive feature according to the classifier. All of the other parameters have less than 10% importance, with the stellar mass ratio and initial eccentricity of the planet being the second-most and third-most important (~8.2%), respectively. τ_a and τ_e are the next two parameters in terms of importance, both with ~5.3%. The planet mass and the ratio of the characteristic damping times are the two least important features with less than 1% importance. The low importance of the time ratio feature is likely due to the fact that the properties that this feature represents in the R_int distribution end up being redundant if the properties represented by τ_a and τ_e have been previously analysed by the classifier.

Considering the high importance of the binary eccentricity with respect to the output MMR of each exoplanet, it is necessary to analyse the distribution of R_int with respect to this parameter. Figure 9 shows the number of simulated exoplanets ejected by each R_int parameterised by e_B. The colour scale shows the percentage representation of each value with respect to the 600 simulations performed for each e_B. From this figure, we note that the higher the binary eccentricity, the higher the value of the predominant capture MMR. Furthermore, it is interesting to note that the label R_int = 4 is mainly obtained for simulations with e_B = {0.01, 0.03, 0.05}, while each remaining R_int has only one principal binary eccentricity associated with it.

Fig. 6 Average values of R obtained for simulations characterised by eB. The y-axis denotes the projection of all the simulations depending on e, q, m, and log10(τa/τe), from top to bottom panels, respectively.

Fig. 7 Confusion matrix of the RFC used, applied to the test data set. The colour bar denotes the percentage that each predicted amount represents normalised to the amount of real values of each label.

Fig. 8 Importance percentage of each feature according to our RFC. The black horizontal lines denote the standard deviation of each estimate.

Fig. 9 Grid with the number of exoplanets according to the binary eccentricity and the final Rint. The colour bar denotes the normalised percentage that each quantity represents out of the total number of simulations for each eccentricity value (600).

5 Discussion

According to the current paradigm, planets do not form in situ around binary stars, but rather migrate from farther out in the disc, where pebble accretion is more favourable. Planets are expected to migrate at a speed proportional to their mass. Remarkably, a slower migration rate makes resonant capture and subsequent ejection more likely.

Thommes & Lissauer (2003) found that two mutually inclined planets migrating in a disc can be captured in resonant inclination when the planetary mass of the outer planet is greater than half that of the inner planet. In our model three-body system, the mass of the planet is negligible in comparison to that of the binary companion. Therefore, we did not consider this resonance configuration in the present work. In principle, it is possible to include an extra term to account for the inclination damping as follows: i(t) = i₀ exp(−t/τ_i), but this requires additional parameterisation. This specific aspect is left for future work.

Martin & Fitzmaurice (2022) showed that while large planets (roughly > 20 M_⊕) may be able to cross the resonances successfully, small planets (roughly < 3 M_⊕) are prone to being captured. These authors considered analytical prescriptions for τ_a and τ_e that modelled Lindblad torques, co-orbital torques, disc surface density, circular gaps in the disc, and stochastic forces to account for the disc turbulence, and demonstrated that the process of resonant ejection of migrating planets may occur in nature. In particular, this mechanism preferentially affects small planets, but is not sufficient to fully explain the dearth of Rp < 3 R_⊕ planets.

Nevertheless, there is likely an important observational bias against the co-existence of small planets due to photometric limits (see discussion by Martin & Triaud 2014).

In this work, we extend the scope of the study by Martin & Fitzmaurice (2022), because their analytical results were developed for larger distances from the binary barycentre (a/a_B ≫ 1). We explored a wider variety of binary eccentricities and damping times. In this sense, our models are compatible with a broader range of circumbinary discs. We find that the 4/1 MMR is the most efficient for resonant capture. However, other resonances can also play a significant role in determining the final parking location of circumbinary planets.

More importantly, circumbinary discs generally form eccentric gaps with a constant precession rate (Thun & Kley 2018; Kley et al. 2019; Penzlin et al. 2021) that may alter the torques from the disc. Surface density (Σ), viscosity (v), disc turbulence (α), and disc metallicity may determine the exoplanet incidence (see e.g. Adibekyan 2019). As shown in Fig. 10, the stellar metallicities of detected CBP systems are statistically lower than the metallicities of single stars harbouring planets. This strongly suggests that we may need to reconsider the typical disc properties around binary stars. To compare with the sample of CBPs, we show only exoplanets around single stars with a radius of less than 12 R_⊕ and periods of less than 300 days. For more realistic discs than the ones considered here, it is likely that the process of MMR capture occurs in any case. Nevertheless, one should not disregard the possibility that, as the planet gains mass and approaches the binary, it switches to type II migration (Armitage 2020). Also, the gas is expected to dissipate after a few million years, which would translate into vanishing disc torques. Lastly, other mechanisms (such as turbulence or mass accretion) might help to avoid the ejection of circumbinary planets.

Fig. 10 Metallicity of stars harbouring known exoplanets compared with the mass of the host star. For binary stars, we assume the total mass (m★ = mB). We identify the groups corresponding to all planets, planets detected by Kepler (including K2), and TOI (including TIC) planets. Exoplanets data obtained from Nasa Exoplanet Archive database.

6 Conclusions

In this work, we analyse the migration of planets with 5, 20, and 93 M_⊕ using Stokes-type forces that mimic planetary migration with a constant rate. Using the mass–radius relation presented in Eq. (1), we estimate that their radii are approximately 1.9, 4.6, and 12.2 R_⊕ (respectively). These values are close to the radius distribution of known CBPs (see Table 1).

As circumbinary disc dissipation timescales are not yet fully understood (Alexander 2012; Shadmehri et al. 2018; Ronco et al. 2021), we considered wide ranges of a-folding and e-folding times, which remain constant throughout the simulation. We also considered a wide range of binary eccentricities. Analysing the final fate of the simulated planets, we find that:

• Circumbinary planets generally migrate along a secular branch until they reach some MMR. At this secular branch, the angle ∆ϖ usually liberates;

• The vast majority of CBPs are ejected at MMRs when migrating towards the binary, with the MMR 4/1 being the one generating the greatest number of ejections;

• The resonant captures always show the resonant angle librating, mainly around 180 degrees $({\psi }_{p/q}^{\left(l\right)}=\left(p+q\right)\lambda -p{\lambda }_B-q\varpi -l\left({\varpi }_B-\varpi \right)$ where q = 1 and l = 1).

Using ensemble learning methods for classification, we draw the following conclusions:

• The binary eccentricity is the most important parameter for determining the (lowest order) MMR at which a CBP may be captured. On the other hand, the planetary mass has low significance;

• It is possible to predict the lowest R_int (defined as the closest integer to the feature R, and associated to the nearest MMR) at which a CBP would be ejected, with a precision of ~81%. To do so, six parameters are required: the binary eccentricity and mass ratio, the planet eccentricity and mass, and the semi-major axis and eccentricity damping times.

Over the last few years, this topic has become an active field of research in planet formation. Indeed, while finishing this study, Fitzmaurice et al. (2022) studied the migration of two circumbinary planets using the approximation of torques induced by a single star. Improvements in our model can be made by assuming some recipes for planetary accretion, which could be important in the runaway regime for large planets (>20 M_⊕) and disc dissipation timescale. The development of more refined models is left for a forthcoming work. In any case, it is of crucial importance to accurately establish how planets migrate within circumbinary discs. This will help us to understand where and for how long circumbinary planets are able to survive around stellar binaries.

Acknowledgements

N-Body computations were performed at Clemente Cluster from IATE, Argentina and at the Mulatona Cluster from the 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. 896319 (SANDS). This research was funded, in part, by ANR (Agence Nationale de la Recherche) of France under contract number ANR-22-ERCS-0002-01. This project has received funding from the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program (grant agreement no. 101042275, project Stellar-MADE).

Appendix A Average ejection times

Fig. A.1 Logarithm of the average ejection times $\left({\log }_{10}\left({\mu }_{t_e{}_j}\right)\right)$ obtained for simulations, characterised by τa and τe. The cells (τa = 106 yr, τe = 105 yr) do not take into account the three non-ejected planets.

Figure A.1 shows a colour scale representing the logarithm of the average ejection times of all simulations in the plane (τ_a, τ_e). For a specific pair of folding times, the values represent the average ejection time and its dispersion within our simulations. We can see that t_ej increases as τ_a and/or τ_e increase(s). In particular, as mentioned in Sect. 2, the mean survival time for each simulation is ~ min(τ_a, τ_e). It is worth mentioning that the errors for this parameter are quite small, giving confidence to the estimates. For example, in (τ_a, τ_e) = (10⁶,10⁶) yr, all simulations, in general, are ejected in $~{10}^{{6.07}_{-0.09}^{+0.07}}$ yr. To corroborate this result, we present Fig. A.2. This figure is analogous to Fig. 6, but shows the averages of t_ej instead of R, and only for the simulations with τ_a = 10⁶ yr and τ_e = 10⁶ yr. All the values are ~t_ej = 10⁶ yr.

Fig. A.2 Logarithm of average ejection times obtained from simulations, characterised by eB. The y-axis denotes the projection of all the simulations depending on e, q, and log10(τa, τe), from top to bottom panels, respectively. In this figure, only the case τa = τe = 106 yr is shown.

Appendix B Resonant angle

Figure B.1 shows the variation of the resonant angle $\left({\psi }_{p/q}^{\left(l\right)}\right)$ and the secular angle (∆ϖ) of the numerical simulations presented in Fig. 3. Here, in the three simulations, we can see the resonant angle circulating until reaching (and being captured by) the corresponding MMR. From this point, the angle starts to librate. By analysing the angle ∆ϖ, we can determine that exoplanets with initial eccentricities e = 0.01 and e = 0.05 migrated while captured by the secular mode (∆ϖ ~ 0). On the other hand, for the extreme case of the exoplanet with initial eccentricity e = 0.1, the semi-major axis and eccentricity dampings are not sufficient for it to be captured by the secular branch. Consequently, in this case, ∆ϖ is continuously circulating.

Fig. B.1 Resonant angle variation ${\psi }_{p/q}^{\left(l\right)}=\left(p+q\right)\lambda -p{\lambda }_B-q\varpi -l\left({\varpi }_B-\varpi \right)$ of the three simulations presented in the top panel of Figure 4. The eccentricities of the planets are e = 0.001 (left), e = 0.01 (centre), and e = 0.1 (right), respectively. Only the three first resonant angles (of each simulation) are shown. The angle ∆ϖ is also shown in the last row.

References

All Tables

All Figures

Fig. 1 Mass–radius relation of transiting circumbinary exoplanets and circumstellar exoplanets. Left panel: radii of exoplanets presented in OC2020 (grey) and CBPs (red/black) as a function of their mass. Their respective error estimations are presented with error bars. In the case of K38, K47b, K413, and K453, we mark their maximum mass reference value with cross symbols (×). The light-blue, orange, and green lines denote the theoretical H2O, Fe, and Earth-like composition lines. The curves outside the box represent the kernel density estimation of the mass (top) and radius (right) of the scattered data points. The data bandwidth for density estimation is calculated with Scott’s Rule (Scott 2015) using the scipy Python package. The vertical (top) and horizontal (right) dashed lines denote the locations of the maximum values in each distribution. All distributions are normalised to their own total counts. Right panel: zoom onto the CBP region shown in the left panel. The exoplanets presented in OC2020 are shown as grey circles, while the CBPs are coloured according to their binary mass. The size of each CBP circle is proportional to the eccentricity of the binary. The red shaded area denotes the mass–radius relation (the red solid line being its reference values) presented in OC2020

In the text

	Fig. 2 Period ratio of known circumbinary planets as a function of their radii. The size of each point represents the planet mass, and the colour scale indicates its binary eccentricity. Horizontal grey lines denote the location of the resonances N/1. The plot is truncated to P/PB|imax = 11.
In the text

Fig. 3 Dynamical maps in the plane (a, e), where each initial condition is integrated for 600 yr and the colour scale corresponds to the ∆e indicator. Solid black lines denote the width of the most prominent N/1 MMR and the dashed vertical lines indicate their nominal location, both calculated using Gallardo et al. (2021). For each of these maps, the binary has q = 0.2 (top frames) and q = 1 (bottom frames), and eB = 0.05 (left frames), eB = 0.3 (centre frames), and eB = 0.5 (right frames). The grey shaded area denotes the stability limit (the blue dot-dashed line being the reference value) calculated by Holman & Wiegert (1999), with its critical eccentricity approximated by Quarles et al. (2018). The solid magenta lines denote the approximation of the ‘capture’ mean eccentricity (Zoppetti et al. 2019). The light-blue dotted line denotes the simulated (a, e) evolution of a planet with m = 5 M⊕ and e = 0.01, setting τa = 105 yr and τe = 106 yr. The red solid line denotes the simple moving average of the planet trajectory, with a window of 20 points (dtw = 105 yr).

In the text

	Fig. 4 Moving average trajectory in the plane (a, e) of three individual simulations with the same τa = 105 yr, τe = 106 yr, q = 0.2, m = 5 M⊕, and eB = 0.05, but different initial eccentricity. Vertical dashed lines denote the location of the 4/1, 5/1, and 6/1 MMRs. The solid magenta lines denote the approximation of the ‘capture mean eccentricity (Zoppetti et al. 2019).
In the text

	Fig. 5 Distribution of R values created with KDE, and a histogram of their closest MMR value (Rint). The values above each Rint bin indicate the fraction of the total 4200 simulations found in that bin as a percentage. The R value KDE distribution is re-scaled for improved visualisation.
In the text

	Fig. 6 Average values of R obtained for simulations characterised by eB. The y-axis denotes the projection of all the simulations depending on e, q, m, and log10(τa/τe), from top to bottom panels, respectively.
In the text

	Fig. 7 Confusion matrix of the RFC used, applied to the test data set. The colour bar denotes the percentage that each predicted amount represents normalised to the amount of real values of each label.
In the text

	Fig. 8 Importance percentage of each feature according to our RFC. The black horizontal lines denote the standard deviation of each estimate.
In the text

	Fig. 9 Grid with the number of exoplanets according to the binary eccentricity and the final Rint. The colour bar denotes the normalised percentage that each quantity represents out of the total number of simulations for each eccentricity value (600).
In the text

	Fig. 10 Metallicity of stars harbouring known exoplanets compared with the mass of the host star. For binary stars, we assume the total mass (m★ = mB). We identify the groups corresponding to all planets, planets detected by Kepler (including K2), and TOI (including TIC) planets. Exoplanets data obtained from Nasa Exoplanet Archive database.
In the text

	Fig. A.1 Logarithm of the average ejection times $\left({\log }_{10}\left({\mu }_{t_e{}_j}\right)\right)$ obtained for simulations, characterised by τa and τe. The cells (τa = 106 yr, τe = 105 yr) do not take into account the three non-ejected planets.
In the text

	Fig. A.2 Logarithm of average ejection times obtained from simulations, characterised by eB. The y-axis denotes the projection of all the simulations depending on e, q, and log10(τa, τe), from top to bottom panels, respectively. In this figure, only the case τa = τe = 106 yr is shown.
In the text

Fig. B.1 Resonant angle variation ${\psi }_{p/q}^{\left(l\right)}=\left(p+q\right)\lambda -p{\lambda }_B-q\varpi -l\left({\varpi }_B-\varpi \right)$ of the three simulations presented in the top panel of Figure 4. The eccentricities of the planets are e = 0.001 (left), e = 0.01 (centre), and e = 0.1 (right), respectively. Only the three first resonant angles (of each simulation) are shown. The angle ∆ϖ is also shown in the last row.

In the text