<figure><img src="image_1.png" /></figure>

HAL
open science

# Circumbinary planets: migration, trapping in mean-motion resonances, and ejection

Emmanuel Gianuzzi, Cristian Giuppone, Nicolás Cuello

► To cite this version:

Emmanuel Gianuzzi, Cristian Giuppone, Nicolás Cuello. Circumbinary planets: migration, trapping in mean-motion resonances, and ejection. Astronomy & Astrophysics - A&A, 2023, 669, pp.A123. 10.1051/0004-6361/202244902 . hal-04970207

HAL Id: hal-04970207

https://hal.science/hal-04970207v1

Submitted on 27 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Circumbinary planets: migration, trapping in mean-motion resonances, and ejection

Emmanuel Gianuzzi<fnref n="1" />,<fnref n="2" />,<fnref n="3" />✉, Cristian Giuppone<fnref n="2" />,<fnref n="3" />✉, and Nicolás Cuello<fnref n="4" />✉

Received 6 September 2022 / Accepted 29 October 2022

## ABSTRACT

**Context.** Most of the planetary systems discovered around binary stars are located at approximately three semi-major axes from the barycentre of their system, curiously close to low-order mean-motion resonances (MMRs). The formation mechanism of these circumbinary planets is not yet fully understood. In situ formation is extremely challenging because of the strong interaction with the binary. One possible explanation is that, after their formation, the interactions between these planets and the surrounding protoplanetary disc cause them to migrate at velocities dependent on the nature of the disc and the mass of the exoplanet. Although extensive data can be obtained with direct hydrodynamical simulations, their computational cost remains too high. On the other hand, the direct *N*-body simulations approach allows us to model a large variety of parameters at much lower cost.

**Aims.** We analyse the planetary migration around a wide variety of binary stars using Stokes-like forces that mimic planetary migration at a constant rate. Our goal is to identify the main parameters responsible for the ejection of planets at different resonances with the inner binary.

**Methods.** We performed 4200 *N*-body simulations with Stokes-like forces and analysed their evolution and outcome as a function of the properties of each system. For each simulated exoplanet, we applied an ensemble learning method for classification in order to clarify the relationship between the inspected parameters and the process of MMR capture.

**Results.** We identify the capture probability for different *N*/1 MMRs, 4/1 being the most prone to capture exoplanets, with 37% probability, followed by MMR 5/1 with ~23% of probability. The eccentricity of the binary is found to be the most important parameter in determining the MMR capture of each circumbinary exoplanet, followed by the mass ratio of the binary and the initial eccentricity of the planet.

**Key words.** binaries: eclipsing – planet-star interactions – methods: numerical – methods: statistical – planets and satellites: dynamical evolution and stability

## 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 \sim 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

<footnote n="1">¹ Facultad de Matemática, Astronomía, Física y Computación (FAMAF), Universidad Nacional de Córdoba (UNC), Córdoba, Argentina
e-mail: egianuzzi@mi.unc.edu.ar</footnote>

<footnote n="2">² Instituto de Astronomía Teórica y Experimental (IATE), CONICET-UNC, Córdoba, Argentina</footnote>

<footnote n="3">³ Observatorio Astronómico de Córdoba (OAC), UNC, Córdoba, Argentina</footnote>

<footnote n="4">⁴ Univ. Grenoble Alpes, CNRS, IPAG/UMR 5274, 38000 Grenoble, France</footnote>
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<fnref n="1" />. 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} \approx 0.06$), which could be a consequence of the interaction with their natal protoplanetary disc. All the planets have radii of greater than $3R_{\oplus}$. Planet masses are typically poorly constrained or compatible with zero. We take metallicities from discovery papers, because we find that the Exoplanet.eu database<fnref n="2" /> has missing values. For instance, the discovery paper of TOI-1338 (Kostov et al. 2020) reports $[Fe/H] = -0.01 \pm 0.05$ dex (measured with HARPS) while its metallicity in the Exoplanet.eu database is $[Fe/H] = -0.4 \pm 0.1$ dex.

<table><caption>**Table 1.** Parameters of the 14 known transiting CBPs.</caption><thead><tr><th>Name</th><th><i>q</i></th><th><i>e</i><sub>bin</sub></th><th><i>R</i><sub>p</sub><br>(<i>R</i><sub>⊕</sub>)</th><th><i>M</i><sub>p</sub><br>(<i>M</i><sub>⊕</sub>)</th><th>Metallicity<br>(dex)</th></tr></thead><tbody><tr><td>K-16</td><td>0.29</td><td>0.16</td><td>8.27</td><td>104.84</td><td>-0.3 [m/H]</td></tr><tr><td>K-34</td><td>0.97</td><td>0.52</td><td>8.38</td><td>69.89</td><td>-0.07 [m/H]</td></tr><tr><td>K-35</td><td>0.91</td><td>0.14</td><td>7.99</td><td>44.78</td><td>-0.34 [m/H]</td></tr><tr><td>K-38</td><td>0.26</td><td>0.10</td><td>4.20</td><td>&lt;122</td><td>-0.11 [m/H]</td></tr><tr><td>K-47b</td><td>0.35</td><td>0.023</td><td>3.05</td><td>&lt;25.77</td><td>-0.25 [m/H]</td></tr><tr><td>K-47d</td><td>0.35</td><td>0.023</td><td>7.04</td><td>19.02</td><td>-0.25 [M/H]</td></tr><tr><td>K-47c</td><td>0.35</td><td>0.023</td><td>4.65</td><td>3.17</td><td>-0.25 [m/H]</td></tr><tr><td>K-64</td><td>0.27</td><td>0.21</td><td>6.10</td><td>168.70</td><td>0.21 [Fe/H]</td></tr><tr><td>K-413</td><td>0.66</td><td>0.037</td><td>4.35</td><td>&lt;168.7</td><td>-0.2 [Fe/H]</td></tr><tr><td>K-453</td><td>0.20</td><td>0.051</td><td>6.30</td><td>&lt;15.88</td><td>0.09 [m/H]</td></tr><tr><td>K-1647</td><td>0.80</td><td>0.16</td><td>11.64</td><td>483.00</td><td>-0.14 [Fe/H]</td></tr><tr><td>K-1661</td><td>0.31</td><td>0.112</td><td>3.87</td><td>17.00</td><td>-0.12 [M/H]</td></tr><tr><td>TOI-1338</td><td>0.29</td><td>0.156</td><td>6.90</td><td>30.20</td><td>0.01 [Fe/H]</td></tr><tr><td>TIC 1729</td><td>0.97</td><td>0.448</td><td>11.25</td><td>942</td><td>0.34 [Fe/H]</td></tr></tbody></table>

**Notes.** K stands for *Kepler*, and TIC 17290098 is identified as TIC-1729. For TIC-1729, we only show the solution named *Family 5* by the authors (Kostov et al. 2021). Discovery papers: 16 (Doyle et al. 2011); 34 & 35 (Welsh et al. 2012); 38 (Orosz et al. 2012a); 47 (Orosz et al. 2012b, 2019); 64 (Schwamb et al. 2013; Kostov et al. 2013); 413 (Kostov et al. 2014); 453 (Welsh et al. 2015); 1647 (Kostov et al. 2016); 1661; (Socia et al. 2020); 1338 (Kostov et al. 2020); 1729 (Kostov et al. 2021).

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_{\oplus}$, with reliable and robust mass and radius measurements. Known CBPs are also included (red markers) for comparison. The distribution of radii of exoplanets is bimodal (David et al. 2021), with peaks at ~2$R_{\oplus}$ and ~13$R_{\oplus}$, while the CBP distribution has a peak close to the valley of this distribution, at ~7$R_{\oplus}$. On the other hand, although the mass distribution of exoplanets can also be categorised as bimodal with peaks at ~9$M_{\oplus}$ and ~260$M_{\oplus}$, 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_{\oplus}$. 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$_2$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) = (0.7 \pm 0.11) \left(\frac{m}{M_{\oplus}}\right)^{(0.63 \pm 0.04)}, \quad (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

<footnote n="1">¹ Obtained from NASA Exoplanet Archive database.</footnote>

<footnote n="2">² http://exoplanet.eu/catalog/</footnote>
<figure><img src="image_3.png" /></figure>

**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 (x). The light-blue, orange, and green lines denote the theoretical H₂O, 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

<figure><img src="image_4.png" /></figure>

**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/P_B|_{max} = 11$.

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 \equiv 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.

## 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_0$ and $m_1$, respectively, and semi-major axis $a_B = 0.1$ AU. The mass of the planet orbiting the binary is equal to $m = 93 M_\oplus$. We choose a Jacobi system for the coordinates. The orbital elements *a*, *e*, *M*, and $\omega$ 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^{-13}$. 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 $\approx 20\,000$ binary periods).

The initial conditions for the circumbinary dynamical maps are coplanar, with mean anomaly $M_B = M = 0^\circ$ and argument of pericentre $\omega_B = \omega = 0^\circ$. To better sample the structure of the resonances, we calculate $\Delta e$ as the amplitude of maximum variation of the orbital eccentricity of the planets during integration, with $\Delta e = (e_{max} - e_{min})$. This indicator correlates with the chaos indicator MEGNO (Cincotta & Simó 2000), and we verified that higher values of $\Delta 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<fnref n="3" />. 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

<footnote n="3">³ www.fisica.edu.uy/~gallardo/atlas/plares.html</footnote>
width of the resonance for coplanar orbits depends on the initial relative value of $\delta\bar{\omega} = \bar{\omega}_B - \bar{\omega} = \omega_B - \omega$. Given that the pericentre of the planet precesses with a short period, the dynamical maps are similar for any initial value of $\delta\bar{\omega}$ 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{d^2\mathbf{r}}{dt^2} = -C(\mathbf{v} - \alpha\mathbf{v}_c), \qquad (2) $$

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

$$ C = \frac{1}{2\tau_a} + \frac{1}{\tau_e}, \qquad \alpha = \frac{2\tau_a}{2\tau_a + \tau_e}, \qquad (3) $$

where $\tau_a$ and $\tau_e$ are the characteristic semi-major axis and eccentricity damping timescales, respectively.

Following Beaugé & 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:

$$ a(t) = a_0 \exp\left(-\frac{t}{\tau_a}\right), \qquad e(t) = e_0 \exp\left(-\frac{t}{\tau_e}\right), \qquad (4) $$

where $a_0$ and $e_0$ are the conditions at the beginning of the integration. The a-folding and e-folding times for $a$ and $e$ are denoted $|\tau_a|$ and $|\tau_e|$, respectively, and have the following expressions:

$$ \tau_a^{-1} = 2C(1 - \alpha), \qquad \tau_e^{-1} = C\alpha. \qquad (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 $\tau_a$, $\tau_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 $\tau_e$ through simulations, some authors use the simple relation

$$ \tau_e = \frac{\tau_a}{K}, \qquad (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 $\tau_a$ and $\tau_e$ independently of each other, and they remain constant throughout all simulations.

<table><caption>**Table 2.** Parameters adopted in the numerical simulations.</caption><thead><tr><th>Parameter</th><th>Adopted values</th></tr></thead><tbody><tr><td><i>q</i></td><td>0.2, 1</td></tr><tr><td><i>e</i><sub>B</sub></td><td>0, 0.01, 0.03, 0.05, 0.1, 0.3, 0.5</td></tr><tr><td><i>e</i></td><td>0, 0.001, 0.01, 0.1</td></tr><tr><td><i>m</i> (<i>M</i><sub>⊕</sub>)</td><td>5, 20, 93</td></tr><tr><td>τ<sub>a</sub> (yr)</td><td>10<sup>3</sup>, 10<sup>4</sup>, 10<sup>5</sup>, 10<sup>6</sup>, 10<sup>7</sup></td></tr><tr><td>τ<sub>e</sub> (yr)</td><td>10<sup>3</sup>, 10<sup>4</sup>, 10<sup>5</sup>, 10<sup>6</sup>, 10<sup>7</sup></td></tr></tbody></table>

**Notes.** Simulations with binary mass ratio $q = m_1/m_0 = 1$ were set up with star masses $m_\star = \{1, 1\} M_\odot$, while the ones with $q = 0.2$ were set up with $m_\star = \{0.833, 0.167\} M_\odot$.

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 ($\tau_a$, $\tau_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(\tau_a, \tau_e)$, we set the output times for each simulation as follows:

$$ dt_{\text{out}} = 0.5 \times 10^{\log_{10}(\min(\tau_e, \tau_a)/\text{yr})-2} \text{yr}. \qquad (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 \times 10^6$ 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 $\Delta 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

$$ a_c = a_B \left( 1.60^{+0.04}_{-0.04} + 5.10^{+0.05}_{-0.05} e_B - 2.22^{+0.11}_{-0.11} e_B^2 \right. \\ \left. + 4.12^{+0.09}_{-0.09} \mu - 4.27^{+0.17}_{-0.17} e_B - 5.09^{+0.11}_{-0.11} \mu^2 \right. \\ \left. + 4.61^{+0.36}_{-0.36} \mu^2 e_B^2 \right), \qquad (8) $$

from Holman & Wiegert (1999), where $\mu = m_1/(m_0 + m_1)$, and using the critical eccentricity approximation,

$$ e_c = 0.8 \left( 1 - \frac{a_c}{a} \right), \qquad (9) $$

from Quarles et al. (2018). These widely used empirical criteria (Yamanaka & Sasaki 2019; Zoppetti et al. 2020; Gallardo et al.
<figure><img src="image_5.png" /><figcaption>**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 e<sub>B</sub> = 0.05 (left frames), e<sub>B</sub> = 0.3 (centre frames), and e<sub>B</sub> = 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<sub>⊕</sub> and e = 0.01, setting τ<sub>a</sub> = 10<sup>5</sup> yr and τ<sub>e</sub> = 10<sup>6</sup> yr. The red solid line denotes the simple moving average of the planet trajectory, with a window of 20 points (dt<sub>w</sub> = 10<sup>5</sup> yr).</figcaption></figure>

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_f^2 + \epsilon^2 \left\langle \frac{T}{L^*} \right\rangle}, \qquad (10) $$

where

$$ e_f = \frac{5a_B e_B (3e_B^2 + 4)(m_0 - m_1)}{8a (3e_B^2 + 2)(m_0 + m_1)} \qquad (11) $$

is the forced eccentricity, and $e^2 \langle T/L^* \rangle$ is an eccentricity term that can be approximated with the following analytical expression (Paardekooper et al. 2012; Zoppetti et al. 2019):

$$ \epsilon^2 \langle T/L^* \rangle = \frac{9}{16} \frac{m_0^2 m_1^2}{(m_0 + m_1)^4} \left(\frac{a_B}{a}\right)^4 \left(1 + \frac{34}{3} e_B^2\right). \qquad (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 \in [0.36, 0.42]$ AU and $e \sim 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 ($\Delta e \le 0.3$) near $a \in [0.35, 0.4]$ AU and $a \in [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_{\oplus}$ with initial eccentricity $e = 0.01$, $\tau_a = 10^5$ yr, and $\tau_e = 10^6$ 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<sup>5</sup> 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(X, N)_j = \frac{\sum_{i=j}^{j+N} X_i}{N}, \qquad (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<sub>B</sub> 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).
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Moving average trajectory in the plane (a, e) of three individual simulations with the same $\tau_a = 10^5$ yr, $\tau_e = 10^6$ yr, $q = 0.2$, $m = 5 M_\oplus$, and $e_B = 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).</figcaption></figure>

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 ($\tau_a = 10^6$ yr, $\tau_e = 10^7$ 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_\oplus$ 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_\oplus, 20 M_\oplus$, and $93 M_\oplus$. 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.

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Distribution of *R* values created with KDE, and a histogram of their closest MMR value (*R*<sub>int</sub>). The values above each *R*<sub>int</sub> 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.</figcaption></figure>

# 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
<table><caption>**Table 3.** Minimum, mean, and maximum index of dispersion of *R* values, calculated for each feature.</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="3">Index of dispersion <i>D</i></th><th rowspan="2">Number of <i>R</i> values</th></tr><tr><th>Min</th><th>Mean</th><th>Max</th></tr></thead><tbody><tr><td><i>e</i><sub>B</sub></td><td>0.044</td><td>0.059</td><td>0.098</td><td>600</td></tr><tr><td><i>q</i></td><td>0.288</td><td>0.328</td><td>0.369</td><td>2100</td></tr><tr><td><i>m</i></td><td>0.325</td><td>0.329</td><td>0.333</td><td>1400</td></tr><tr><td><i>e</i></td><td>0.307</td><td>0.321</td><td>0.336</td><td>1050</td></tr><tr><td><i>τ</i><sub>a</sub></td><td>0.275</td><td>0.325</td><td>0.368</td><td>840</td></tr><tr><td><i>τ</i><sub>e</sub></td><td>0.268</td><td>0.322</td><td>0.359</td><td>840</td></tr><tr><td><i>τ</i><sub>a</sub>/<i>τ</i><sub>e</sub></td><td>0.273</td><td>0.317</td><td>0.340</td><td>168</td></tr></tbody></table>

**Notes.** The last column shows the number of simulations used to calculate each index of dispersion *D*.

of each bin indicate the percentage of ejected planets for a given *R*<sub>int</sub>. 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<fnref n="4" /> 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. *τ*<sub>a</sub>/<i>τ</i><sub>e</sub> is also considered as an additional (derived) parameter. The distribution of *R* values for each *e*<sub>B</sub> value has the lowest dispersion index of all the features. Therefore, *e*<sub>B</sub> 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<sub>10</sub>(*τ*<sub>a</sub>/<i>τ</i><sub>e</sub>) (bottom frame) on the average *R* value (*μ*<sub>R</sub> = <i>R̄</i>) and its dispersion (*σ*<sub>R</sub>). All grids have *e*<sub>B</sub> 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 *μ*<sub>R</sub> between planets with low eccentricity (*e* < 0.1).

From the second panel, it is possible to infer that the dependence of *μ*<sub>R</sub> on *q* is weakly related to the binary eccentricity. For a large value of the mass ratio, the increase rate of *μ*<sub>R</sub> with respect to *e*<sub>B</sub> 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) *μ*<sub>R</sub>/<i>e</i><sub>B</sub>, while that for simulations with *q* = 0.2 is ~ (6.49 ± 0.60) *μ*<sub>R</sub>/<i>e</i><sub>B</sub>.

From the third panel, the value of *μ*<sub>R</sub> 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<sub>10</sub>(*τ*<sub>a</sub>/<i>τ</i><sub>e</sub>) and *μ*<sub>R</sub>. The value of *μ*<sub>R</sub> decreases as the value of |log<sub>10</sub>(*τ*<sub>a</sub>/<i>τ</i><sub>e</sub>)| increases, regardless of whether *τ*<sub>a</sub> is higher or lower than *τ*<sub>e</sub>. This effect is more noticeable as the binary eccentricity increases. It is interesting to note that, although the values of *μ*<sub>R</sub> are similar, their dispersion values *σ*<sub>R</sub> are slightly higher for the simulation groups with *τ*<sub>a</sub> ≤ *τ*<sub>e</sub>.

## 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 (Ivezić 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`<fnref n="5" /> 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

<footnote n="4">⁴ In this section we prefer to use the term ‘feature’ instead of ‘parameter’, as this is typical machine learning terminology.</footnote>

<footnote n="5">⁵ https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 6.** Average values of *R* obtained for simulations characterised by *e*<sub>B</sub>. The *y*-axis denotes the projection of all the simulations depending on *e*, *q*, *m*, and log<sub>10</sub>(&tau;<sub>a</sub>/&tau;<sub>e</sub>), from top to bottom panels, respectively.</figcaption></figure>

each set are [$e_B$, $e$, $\tau_a$, $\tau_e$, $q$, $m$, $\tau_a/\tau_e$], and their target is the integer closest to each value $R$, named $R_{int} \equiv \text{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

<figure><img src="image_9.png" /><figcaption>**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.</figcaption></figure>

<table><caption>**Table 4.** Classification report of the RFC used.</caption><thead><tr><th>R<sub>int</sub></th><th>Precision</th><th>Recall</th><th>F1-score</th><th>Support</th></tr></thead><tbody><tr><td>3</td><td>0.91</td><td>0.82</td><td>0.86</td><td>98</td></tr><tr><td>4</td><td>0.82</td><td>0.89</td><td>0.85</td><td>402</td></tr><tr><td>5</td><td>0.75</td><td>0.71</td><td>0.73</td><td>235</td></tr><tr><td>6</td><td>0.78</td><td>0.67</td><td>0.72</td><td>145</td></tr><tr><td>7</td><td>0.69</td><td>0.75</td><td>0.72</td><td>121</td></tr><tr><td>8</td><td>0.72</td><td>0.72</td><td>0.72</td><td>46</td></tr><tr><td>9</td><td>1.00</td><td>0.67</td><td>0.80</td><td>3</td></tr><tr><td>Macro average</td><td>0.81</td><td>0.75</td><td>0.77</td><td>1050</td></tr><tr><td>Weighted average</td><td>0.79</td><td>0.79</td><td>0.79</td><td>1050</td></tr></tbody></table>

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
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Importance percentage of each feature according to our RFC. The black horizontal lines denote the standard deviation of each estimate.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Grid with the number of exoplanets according to the binary eccentricity and the final $R_{int}$. The colour bar denotes the normalised percentage that each quantity represents out of the total number of simulations for each eccentricity value (600).</figcaption></figure>

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. $\tau_a$ and $\tau_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 $\tau_a$ and $\tau_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.

## 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_0 \exp(-t/\tau_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 $\tau_a$ and $\tau_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 $R_p < 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 \gg 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 ($\Sigma$), viscosity ($\nu$), disc turbulence ($\alpha$), 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
<figure><img src="image_12.png" /><figcaption>**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_{\star} = m_{\odot}$). 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.</figcaption></figure>

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_{\oplus}$ 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.

## 6. Conclusions

In this work, we analyse the migration of planets with 5, 20, and 93 $M_{\oplus}$ 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_{\oplus}$ (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 $\Delta\sigma$ 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 librat- ing, mainly around 180 degrees ($\psi_{p/q}^{(l)} = (p+q)\lambda - p\lambda_{\text{B}} - q\sigma - l(\omega_{\text{B}} - \sigma)$ 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_{\text{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_{\oplus}$) 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).

## References

Adibekyan, V. 2019, *Geosciences*, **9**, 105

Alexander, R. 2012, ApJ, **757**, L29

Armitage, P. J. 2020, *Astrophysics of Planet Formation*, 2nd edn. (Cambridge University Press, UK)

Audenaert, J., Kuszlewicz, J. S., Handberg, R., et al. 2021, AJ, **162**, 209

Beauge, C., & Ferraz-Mello, S. 1993, Icarus, **103**, 301

Beaugé, C., Michtchenko, T. A., & Ferraz-Mello, S. 2006, MNRAS, **365**, 1160

Bengfort, B., & Bilbro, R. 2019, J. Open Source Softw., **4**

Breiman, L. 2001, Mach. Learn., **45**, 5

Bromley, B. C., & Kenyon, S. J. 2015, ApJ, **806**, 98

Chrenko, O., Brož, M., & Nesvorný, D. 2018, ApJ, **868**, 145

Cincotta, P. M., & Simó, C. 2000, A&AS, **147**, 205

Cox, D. R., & Lewis, P. A. W. 1966, *The Statistical Analysis of Series of Events* (Springer)

Cresswell, P., & Nelson, R. P. 2008, A&A, **482**, 677

David, T. J., Contardo, G., Sandoval, A., et al. 2021, AJ, **161**, 265

Doyle, L. R., Carter, J. A., Fabrycky, D. C., et al. 2011, Science, **333**, 1602

Esmer, E. M., Baştürk, Ö., Selam, S. O., & Aliş, S. 2022, MNRAS, **511**, 5207

Fitzmaurice, E., Martin, D. V., & Fabrycky, D. C. 2022, MNRAS, **512**, 5023

Gallardo, T., Beaugé, C., & Giuppone, C. A. 2021, A&A, **646**, A148

Giuppone, C. A., Rodríguez, A., Michtchenko, T. A., & de Almeida, A. A. 2022, A&A, **658**, A99

Hirsh, K., Price, D. J., Gonzalez, J.-F., Ubeira-Gabellini, M. G., & Ragusa, E. 2020, MNRAS, **498**, 2936

Holman, M. J., & Wiegert, P. A. 1999, AJ, **117**, 621

Ivezić, Ž., Connelly, A. J., VanderPlas, J. T., & Gray, A. 2014, Statistics, Data Mining, and Machine Learning in Astronomy, Princeton Series in Modern Observational Astronomy (Princeton University Press)

Kley, W., & Haghighipour, N. 2014, A&A, **564**, A72

Kley, W., & Haghighipour, N. 2015, A&A, **581**, A20

Kley, W., & Nelson, R. P. 2012, ARA&A, **50**, 211

Kley, W., Papaloizou, J. C. B., & Ogilvie, G. I. 2008, A&A, **487**, 671

Kley, W., Thun, D., & Penzlin, A. B. T. 2019, A&A, **627**, A91

Kostov, V. B., McCullough, P. R., Hinse, T. C., et al. 2013, ApJ, **770**, 52

Kostov, V. B., McCullough, P. R., Carter, J. A., et al. 2014, ApJ, **784**, 14

Kostov, V. B., Orosz, J. A., Welsh, W. F., et al. 2016, ApJ, **827**, 86

Kostov, V. B., Orosz, J. A., Feinstein, A. D., et al. 2020, AJ, **159**, 253
Kostov, V. B., Powell, B. P., Orosz, J. A., et al. 2021, AJ, 162, 234

Lee, M. H., & Peale, S. J. 2002, ApJ, 567, 596

Lu, Y., Angus, R., Agüeros, M. A., et al. 2020, AJ, 160, 168

Martin, D. V. 2018, Populations of Planets in Multiple Star Systems (Springer), 156

Martin, D. V., & Fitzmaurice, E. 2022, MNRAS, 512, 602

Martin, D. V., & Triaud, A. H. M. J. 2014, A&A, 570, A91

Martin, D. V., Triaud, A. H. M. J., Udry, S., et al. 2019, A&A, 624, A68

Meschiari, S. 2012, ApJ, 761, L7

Moriwaki, K., & Nakagawa, Y. 2004, ApJ, 609, 1065

Mutter, M. M., Pierens, A., & Nelson, R. P. 2017, MNRAS, 469, 4504

Nelson, R. P. 2003, MNRAS, 345, 233

Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012a, ApJ, 758, 87

Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012b, Science, 337, 1511

Orosz, J. A., Welsh, W. F., Haghighipour, N., et al. 2019, AJ, 157, 174

Otegi, J. F., Bouchy, F., & Helled, R. 2020, A&A, 634, A43

Paardekooper, S.-J., Leinhardt, Z. M., Thébault, P., & Baruteau, C. 2012, ApJ, 754, L16

Paardekooper, S.-J., Dong, R., Duffell, P., et al. 2022, ArXiv e-prints [arXiv:2203.09595]

Penzlin, A. B. T., Ataiee, S., & Kley, W. 2019, A&A, 630, A1

Penzlin, A. B. T., Kley, W., & Nelson, R. P. 2021, A&A 645, A68

Pierens, A., & Nelson, R. P. 2008, A&A, 478, 939

Pierens, A., & Nelson, R. P. 2013, A&A, 556, A134

Pierens, A., McNally, C. P., & Nelson, R. P. 2020, MNRAS, 496, 2849

Poblete, P. P., Cuello, N., & Cuadra, J. 2019, MNRAS, 489, 2204

Quarles, B., Satyal, S., Kostov, V., Kaib, N., & Haghighipour, N. 2018, ApJ, 856, 150

Ragusa, E., Dipierro, G., Lodato, G., Laibe, G., & Price, D. J. 2017, MNRAS, 464, 1449

Ragusa, E., Alexander, R., Calcino, J., Hirsh, K., & Price, D. J. 2020, MNRAS, 499, 3362

Rein, H. 2012, MNRAS, 427, L21

Ribas, Á., Merín, B., Bouy, H., & Maud, L. T. 2014, A&A, 561, A54

Ronco, M. P., Schreiber, M. R., Giuppone, C. A., et al. 2020, ApJ, 898, L23

Ronco, M. P., Guilera, O. M., Cuadra, J., et al. 2021, ApJ, 916, 113

Schlecker, M., Pham, D., Burn, R., et al. 2021, A&A, 656, A73

Schwamb, M. E., Orosz, J. A., Carter, J. A., et al. 2013, ApJ, 768, 127

Schwarz, R., Funk, B., Zechner, R., & Bazsó, Á. 2016, MNRAS, 460, 3598

Scott, D. W. 2015, Multivariate Density Estimation: Theory, Practice, and Visualization (Hoboken, NJ: John Wiley and Sons)

Secunda, A., Bellovary, J., Mac Low, M.-M., et al. 2019, ApJ, 878, 85

Secunda, A., Bellovary, J., Mac Low, M.-M., et al. 2020, ApJ, 903, 133

Shadmehri, M., Ghoreyshi, S. M., & Alipour, N. 2018, ApJ, 867, 41

Silsbee, K., & Rafikov, R. R. 2015, ApJ, 798, 71

Socia, Q. J., Welsh, W. F., Orosz, J. A., et al. 2020, AJ, 159, 94

Sutherland, A. P., & Kratter, K. M. 2019, MNRAS, 487, 3288

Thommes, E. W., & Lissauer, J. J. 2003, ApJ, 597, 566

Thun, D., & Kley, W. 2018, A&A, 616, A47

Thun, D., Kley, W., & Picogna, G. 2017, A&A, 604, A102

Ulmer-Moll, S., Santos, N. C., Figueira, P., Brinchmann, J., & Faria, J. P. 2019, A&A, 630, A135

Welsh, W. F., Orosz, J. A., Carter, J. A., et al. 2012, Nature, 481, 475

Welsh, W. F., Orosz, J. A., Short, D. R., et al. 2015, ApJ, 809, 26

Yamanaka, A., & Sasaki, T. 2019, Earth Planets Space, 71, 82

Zoppetti, F. A., Beaugé, C., & Leiva, A. M. 2018, MNRAS, 477, 5301

Zoppetti, F. A., Beaugé, C., & Leiva, A. M. 2019, in Journal of Physics Conference Series, 1365, 012029

Zoppetti, F. A., Leiva, A. M., & Beaugé, C. 2020, A&A, 634, A12
## Appendix A: Average ejection times

<figure><img src="image_13.png" /><figcaption>**Fig. A.1.** Logarithm of the average ejection times ($\log_{10}(\mu_{t_{e_j}})$) obtained for simulations, characterised by $\tau_a$ and $\tau_e$. The cells ($\tau_a = 10^6$ yr, $\tau_e = 10^5$ yr) do not take into account the three non-ejected planets.</figcaption></figure>

Figure A.1 shows a colour scale representing the logarithm of the average ejection times of all simulations in the plane ($\tau_a, \tau_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_{e_j}$ increases as $\tau_a$ and/or $\tau_e$ increase(s). In particular, as mentioned in Sect. 2, the mean survival time for each simulation is ~min($\tau_a, \tau_e$). It is worth mentioning that the errors for this parameter are quite small, giving confidence to the estimates. For example, in ($\tau_a, \tau_e$) = ($10^6, 10^6$) 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_{e_j}$ instead of R, and only for the simulations with $\tau_a = 10^6$ yr and $\tau_e = 10^6$ yr. All the values are ~$t_{e_j} = 10^6$ yr.

<figure><img src="image_14.png" /><figcaption>**Fig. A.2.** Logarithm of average ejection times obtained from simulations, characterised by $e_B$. The y-axis denotes the projection of all the simulations depending on $e, q$, and $\log_{10}(\tau_a/\tau_e)$, from top to bottom panels, respectively. In this figure, only the case $\tau_a = \tau_e = 10^6$ yr is shown.</figcaption></figure>

## Appendix B: Resonant angle

Figure B.1 shows the variation of the resonant angle ($\psi_{p/q}^{(l)}$) and the secular angle ($\Delta\omega$) 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 $\Delta\omega$, we can determine that exoplanets with initial eccentricities $e = 0.01$ and $e = 0.05$ migrated while captured by the secular mode ($\Delta\omega \sim 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, $\Delta\omega$ is continuously circulating.
<figure><img src="image_15.png" /><figcaption>**Fig. B.1.** Resonant angle variation $\psi_{p/q}^{(l)} = (p+q)\lambda - p\lambda_B - q\omega - l(\omega_B - \omega)$ 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 $\Delta\omega$ is also shown in the last row.</figcaption></figure>