© The Authors 2023

1 Introduction

Until recently, coorbital configurations or 1:1 resonances could only be detected in small objects in the Solar System (e.g. Janus and Epimetheus) possibly also due to detection biases (Giuppone et al. 2012). Theoretically, around a single star, Leleu et al. (2018) found stable 1:1 resonances for planets with masses of 10⁻⁵ M_⋆ for orbit eccentricities up to 0.5. Raymond et al. (2023) recently demonstrated that even multi-planet coorbital configurations can be stable for low masses. Cresswell & Nelson (2006) found in single-star systems that the migration of protoplanets can lead to 1:1 resonant horseshoe orbits. However, to achieve stability in the single star case, the resonance with a larger additional planet is needed (Cresswell & Nelson 2009; Leleu et al. 2019). Without the additional resonant planet continued migration leads to increasing libration and breaking the resonance. Recently, Balsalobre-Ruza et al. (2023) found evidence for accumulated material at the Lagrange L5 point of the orbit of the embedded planet PDS 70 b. This could indicate the formation of coor-bital Trojan objects in protoplanetary systems as theoretically predicted.

In binary star systems, migration of the planet through the circumbinary disc can be halted by the eccentric inner cavity in the disc caused by the binary (Mutter et al. 2017; Kley et al. 2019; Penzlin et al. 2021). The eccentricity of the inner disc near the cavity refills planet-opened gaps and thereby hinders planets from reaching the inner cavity. The inner circumbinary disc is most eccentric for low (<0.05) and high (>0.3) binary eccentricities (Thun et al. 2017; Ragusa et al. 2020). This leads to the condition that allows a 1:1 resonance to form in a circumbinary disc of nearly circular binaries. In Penzlin et al. (2019), we found that a 1:1 resonance can form between comparable mass planets of ≈10⁻⁵−10⁻⁴ M_bin inside the circumbinary disc. The migration of the planets is halted at comparable positions due to the structure of the inner disc and the damping forces of the eccentric disc allowing them to enter a coorbital configuration. The 1:1 resonance in the eccentric orbit, which resulted from the shape of the eccentric disc, remains stable even without the disc in n-body simulations. The n-body models, including a prescribed disc by Coleman et al. (2023), also find instances of dynamically created coorbital planet configurations for Kepler-16 and Kepler-34 analogues (Doyle et al. 2011; Welsh et al. 2012). During the embedded phase, the structure of the circumbinary disc leads to slow planet migration, exposing them to potentially unstable resonances with the binary on their orbit as studied in Martin & Fitzmaurice (2022) for single planets.

In the case of a coorbital resonance, the destabilizing resonances may be even more important as the two planets librate around their common orbit and thereby enter mean motion resonant orbits with the binary even at distances farther from the resonant location (see also Sutherland & Kratter 2019, for resonant multi-planets). Additionally, the resonant location widens with increased orbit eccentricity (Mudryk & Wu 2006). This prompts the question of which other circumbinary orbits are theoretically stable for 1:1 resonant planets.

In the last decade, a dozen close circumbinary planets have been detected with the latest confirmed observation of TIC 172900988b (Kostov et al. 2021) and TOI-1338 b&c (Standing et al. 2023). Most detected single circumbinary planets show comparable orbits in units of binary distance close to the instability limit (Dvorak 1986; Holman & Wiegert 1999). Studies by Popova & Shevchenko (2016); Quarles et al. (2018) and Chen et al. (2019) have shown that also eccentric single planets close to the binary can be stable, and they described how the closest stable orbit depends on the planet’s eccentricity.

One special case is the Kepler-47 system that hosts three circumbinary planets, demonstrating that multi-planet configurations around binaries are possible (Orosz et al. 2019). In Penzlin et al. (2019), we found stable 1:1 resonant configurations for Kepler-47 and -413 analogues in hydrodynamic and n-body simulations. With this idea, Sudol & Haghighipour (2021) constrained the observability of resonant planets.

A theoretical study by Michalodimitrakis & Grigorelis (1989) found only stable 1:1 resonances around binaries for zero-mass planets. However, since we already identified at least one stable orbit around a circular binary system with planet masses of ≤0.21 M_jup, for this study we investigated stable planet orbits over the full range of planetary eccentricities with rebound n-body simulations (Rein & Liu 2012; Rein & Spiegel 2015) to identify the stability limit of 1:1-resonant, circumbinary planets with a Kepler-413-like host binary and for a range of binary parameters relevant to the known circumbinary systems.

Section 2 explains the n-body model used to analyse the orbit stability. In Sect. 3, Kepler-413 is used as an example to test the best set-up, explore the dynamics of the planet systems, and establish a simple first order stability fitting function. We discuss the resulting parameter range of 1:1 resonant stable planet orbits for different binaries in Sect. 4. Finally, we discuss the results in Sect. 5 and summarize the findings in Sect. 6.

Fig. 1 Initial configuration of the stars (S1 and S2) and planets (P1 and P2) with the relative argument of periastron Δω between planet and binary orbits.

2 Models

In this study, we investigate a binary system hosting circumbi-nary planets in an initial coorbital configuration (see Fig. 1). We compare these to systems containing only a single circumbinary planet. The initial orbit parameters of the n-body problem, the variations of parameters, and the simulation setup are defined as follows.

2.1 Binary parameters

We use the Kepler-413 binary as a fiducial case to go onto exploring a range of binary parameters. We base the fiducial simulation on the parameters found by Kostov et al. (2014), therefore the binary mass ratio is µ_bin = M₂/M_bin = 0.4 using the primary mass M₁ and total binary mass M_bin and the binary eccentricity is e_bin = 0.04. Mass and distance are scalable with the binary parameters used, so we use the binary semi-major axis 1 a_bin for one standard length unit and 1 M_bin for the standard mass unit. The gravitational constant is normalized to 1, which results in a binary period of T_bin = 2π, which we use as the standard time unit. The eccentricity of the planets varies between 0.0–0.9 in the stability analysis. The setup is summarized in Table 1 as fiducial.

Investigating the binary parameters, the mass ratio µ_bin takes the values [0.1, 0.2, 0.3, 0.4, 0.5] and the eccentricity e_bin [0.0, 0.125, 0.25, 0.375, 0.5] at the same time. This creates a 5 × 5 parameter space to understand the rudimentary effect of the binary on the planet’s stable orbits (see ‘par’ in Table 1). The binary orbit starts at its apocentre position.

2.2 Planet parameters

The planetary orbits are initialized relative to the binary centre of mass so that the relative argument of pericentre Δω = ω_bin − ω_p defines the relative alignment of the two orbits (see Fig. 1). The argument of pericentre and the orbital eccentricity of both coor-bital planets are identical, where the first planet starts at the pericentre of the common orbit and the second one at the apoc-entre. Most of the time this results in planets on a horseshoe orbit configuration. Other simulations showed that different values for the initial true anomalies of the planets may result in a different kind of coorbital configuration (Tadpole orbit). However, studying our fiducial system revealed the horseshoe configuration to be the most interesting, therefore we focussed the study on that (see also Sect. 5).

Since the potential of the binary differs from a single object, the initial velocities and semi-major axes of the planets are constructed considering the centre of mass rather than the exact potential and, therefore, can differ from the dynamical ideal solution. However, the difference depends on Δω and is minor. We discuss the influence in the two planets case in Sect. 3. To have a conservative estimate, we generally choose the least stable initial angle of π/2 according to our parameter studies for planet eccentricities of 0.1, 0.3, and 0.5, which is also in agreement with Quarles et al. (2018).

As an example, the e_p = 0.3 coorbital case is displayed in Fig. 2. The effect is more significant in the coorbital configuration.

The planet masses are equal and have a combined mass of 10⁻⁴ M_bin, or ~30 M_⊕. This mass is typical considering the planets observed in a binary system from the Kepler mission. To resolve the stable orbits we simulated all orbits with a spacing of Δr_peri = 0.01 a_bin starting at a clearly unstable orbit ≥2 a_bin.

Fig. 2 Stability of coorbital planets with an orbital eccentricity of 0.3 depending on the initial angle between binary periapsis and planet-orbit periapsis. White colour shows stable simulations and the colours indicate the survival time of the coorbital configuration. The dashed and dotted lines mark the 1:x and 2:x planet-binary resonances respectively.

2.3 Simulation setup

To simulate each binary system, we use the rebound n-body code by Rein & Liu (2012) with the IAS15-integrator (Rein & Spiegel 2015) for a 2D setup¹. To evaluate the stability, we define the following termination conditions: The planetary semi-major axis changes substantially. For a single planet, the maximum allowed variation is 20% of its initial value. For 1:1 planets, the maximum allowed difference between planetary semi-major axes is 20% of the inner planet’s semi-major axis evaluated at each time step.

Coorbital configurations must maintain a nearly identical semi-major axis as a condition for the 1:1 resonance. As objects on a horseshoe orbit can get closer than 0.1 a_bin despite the orbit staying stable, we do not use a close encounter termination condition. However, as all objects are considered point masses, any collision would lead to an ejection, violating the above condition. We consider a simulation stable if it reaches 10⁵ T_bin without termination.

3 Stability of single and 1:1 resonant planets around Kepler-413

First, we highlight the stability of single and 1:1 resonant planet systems around a binary similar to Kepler-413, because such a low-eccentric binary causes an eccentric disc (Kley et al. 2019; Ragusa et al. 2020; Muñoz & Lithwick 2020; Penzlin et al. 2021) that more easily traps planets into the 1:1 resonance (Penzlin et al. 2019). With this example, we explore some of the caveats and the properties to consider for the dynamics of the planet pair to choose a setup that does not overestimate the stability.

For the single planet configuration, the planet is added in a Keplerian orbit around the centre of mass of the binary. To initialize the complete 1:1 resonance configuration, we add the second planet in the same Keplerian orbit calculated around the centre of mass of the entire single planet configuration, but with a true anomaly of π, meaning in its apocentre. The resulting initial coordinates for the example of a_p,ini = 5.45 are given in Table 2.

Due to the sequential initialization of the planets, the mean orbit can vary in semi-major axis from the initially set value axis by up to 2% as shown in Fig. 3. The alignment becomes a factor due to the non-perfect initial configuration of the planets that also does not include the higher order structure of the binary potential. As discussed before, the initial relative argument of periastron between eccentric planet orbit and eccentric binary orbit affects stability as shown in the example case in Fig. 2 for e_p = 0.3. Compared to a single planet, the coorbital configuration is much less stable and more likely to be interrupted by binary-planet resonances, since the libration of the 1:1 resonant planets leads to natural variations in the semi-major axis of both planets, as described in Mudryk & Wu ( 006). In Fig. 2 the lines mark the positions of the 1:x and 2:x resonances with the binary. The inner edge of some stability regions (light colours in Fig. 2) coincides with 1:x resonances, while the 2:x resonances limit the outer edges. The width of these resonant instability regions decreases for increasing semi-major axis as the resonances also become narrow, weakened by the overlap of more higher order resonances and the reduced gravitational potential of the fast-moving stars on the planets. Thus orbits further out become stable.

In general, the stability variations due to the misalignment between binary and planet orbits Δω is naturally decreased by reducing the eccentricity of the planet’s orbit. Varying the initial orbital distance leads to varying the orbital stability. For less eccentric planet orbits, the closest distance between a planet and binary varies less with Δω and, in the case of a circular planet orbit, the closest distance remains constant at the difference between the planet orbit radius and binary pericentre. We confirm the findings by Quarles et al. (2018) that initially aligned configurations allow more systems to settle into apparent stable states than those with perpendicular oriented arguments of periastron. So for Δω = π/2 the unstable orbits are maximal, therefore we use Δω = π/2 in the eccentricity runs to have a conservative measure of the stability of the orbits in the parameter study.

The relative orbits between the planets are generally horseshoe orbits as shown in Fig. 4 for the example of e_p = 0.3, a_p = 5.7 a_bin. The planets librate around the Lagrange point locations L4 and L5 and exchange their angular momentum with the orbit swap (see Fig. 3). The pattern is more complex due to the effects of the binary; however, it is periodic since the colour scale from yellow to blue tracks the linear time progression of 10⁵ T_bin. To investigate the dynamics of the planets, we compare our findings with the expected precession period from the restricted three-body problem. Figure 5 shows the precession of the common planet orbit for a planet eccentricity of e_p = 0.3. The time is normalized to the theoretical precession of the restricted three-body problem derived from Moriwaki & Nakagawa (2004) $$T_{\text{prec }}=\frac{4}{3}\frac{{(q+1)}^2}{q}{a}^{3.5}\frac{{\left(1-e_{\text{p}}^2\right)}^2}{1+1.5e_{\text{bin }}^2}{T}_{\text{bin }}.$$(1)

In the selected cases, the precession of the resonant planets is systematically slightly slower by ≈50–70 T_bin for stable orbits. However, this precession period is still close to the prediction of the restricted three-body problem. This does not change by reducing the planet masses to 10⁻⁶ M_bin.

To compare to previous work we simulate a range of single circumbinary planet systems and vary the initial planetary eccentricity and semi-major axis. Figure 6 shows the instability transition with increasing planetary eccentricity and correspondingly increasing pericentre distance. The non-blue region denotes stable initial conditions for a single planet (up to 10⁵ T_bin). Figure 6 demonstrates that with increasing eccentricity the first stable orbit moves further out for the single planet. By plotting the pericentre distance of the first stable planet orbit a linear trend can be found, with the exception of very small planet eccentricities. Due to the 1:4 resonance, the values for e_p < 0.2 do not follow the overall upward trend of instability. The 1:4 resonance corresponds to the lowest stable region in Fig. 6 surrounding the lowest black dashed line plotted. The effect of the 1:x resonances in the single planet case beyond 1:8 becomes hardly noticeable. The critical semi-major axis for the circular planet of 2.3 a_bin compares well with the value for a circular planet around a circular binary in Holman & Wiegert (1999). Figure 7 shows the fitted trend of stability for the single planet in blue. Due to resonance some orbits further out are unstable (see Sect. 4), however, for the more general trend using the first stable orbit is sufficient. The fitted stability limit is $$r_{\text{c},\text{\hspace{0.17em}peri }}=\left[(2.08\pm 0.10)e_{\text{p}}+(1.91\pm 0.05)\right]a_{\text{bin}},$$(2)

using standard deviation errors.

The simple linear fit describes the stability behaviour of a single planet around a similar to Kepler-413 fairly well for e_p > 0.2. However, the fit ignores the resonant unstable lines. Especially the 1:4 and 1:3 resonant instability decrease the stability for smaller planet eccentricities (e_bin < 0.1) below a linear fit. The fit also slightly overestimates the stability for the most eccentric systems, which is realized in other works making use of a quadratic term (e.g. Mardling 2001).

Now focussing on the results of the coorbital planet configuration, Fig. 6 shows the stability of orbits for different planet eccentricities and pericentre positions up to an eccentricity of 0.9 in purple. Stable orbits exist for all eccentricities if the orbit is sufficiently large.

For the 1:1 resonant planets as well, the stable orbits follow a linear trend, plotted in Fig. 7. The linear fitting function is given by: $$r_{\text{c},\text{\hspace{0.17em}peri }}=\left[(2.90\pm 0.15)e_{\text{p}}+(2.46\pm 0.08)\right]a_{\text{bin}}.$$(3)

The offset of this function here is larger than for the single planet, meaning that stable orbits in the case of 1:1 resonant planets around a binary start >0.5 a_bin further out, and the function for 1:1 resonant planets is steeper with planet eccentricity. So single planets with eccentric orbits can survive much closer to the binary than 1:1 resonant planet configurations. This is a result of the broadened libration of the planet orbits due to the added coorbital angular momentum exchange.

However, stable 1:1 resonant planets around binaries are possible within a pericentre of r_peri < 5.6 a_bin. This means that the planets still survive on orbits that pass by close to the binary even for extremely eccentric 1:1 configurations.

Fig. 3 Semi-major axis evolution of the two planets for a time window of 4000 binary orbits after initialization, showing the libration of the planets around a common orbit. The initial semi-major axis was 5.45 abin. However, the average semi-major axis (black dashed line) during the simulation denoted in the top right corner is slightly smaller.

Fig. 4 Plot of the relative mean longitude and semi-major axis between the coorbital planets at ap = 5.7 abin with ep = 0.3, indicating that the planets’ relative motion describes a stable horseshoe orbit around the L4 and L5 points marked with red circles, without orbit crossing. The colour progresses over 105 Tbin from yellow to dark blue.

Fig. 5 Relative argument of periastron of the common planets orbit and the binary over the precession rate normalized to the theoretical precession rate from the restricted three-body problem.

Fig. 6 Map of the survival of single and double (coorbital) planet systems with varying planet eccentricities ep. Blue is unstable for all, purple is only stable for single planets but not the coorbital configuration. The black lines mark l:x resonances between x = 4−13 from shallowest to steepest.

Fig. 7 Linear fit of the pericentre of the first stable orbits in Fig. 6. The values close to ep = 0 diverge from the plot due to a unstable region around the 1:4 resonance.

4 Range of stable 1:1 resonance

After using this specific example of a Kepler-413-like system, we vary the binary parameters to compare the results to the general stability of planets around binaries as investigated before for the single planet case (Dvorak 1986; Holman & Wiegert 1999; Chavez et al. 2015; Popova & Shevchenko 2016; Quarles et al. 2018) and create a simple model of the stability limits around binary stars. For this, we change the mass fraction µ = {0.1, 0.2, 0.3, 0.4, 0.5}, while also changing the binary eccentricity e_bin = {0, 0.125, 0.25, 0.327, 0.5}. We limited the range of binary eccentricities in this study to the range of the observed host binary stars, with Kepler-34 being the most eccentric, planet-hosting binary with e_bin = 0.52 (Welsh et al. 2012). In Fig. 8 we show the results of the orbital stability for the whole data set. The colour marks the logarithmic survival time of the system, where white systems survived 10⁵ T_bin. More equal mass binaries reduce the survival of more eccentric planet configurations. More eccentric binary stars shift the first stable orbit in general further out as is also suggested in other work on single planets (e.g. Chavez et al. 2015; Popova & Shevchenko 2016; Quarles et al. 2018). Interestingly, the moderately binary eccentricities lead to unstable regions further out after the first stable orbit for configurations with planets more eccentric than 0.5. For eccentric binary systems, the resonant orbits cause even more pronounced unstable feature lines.

Comparing the single planet and coorbital planet configurations, the slope of the unstable region is steeper for coorbital configurations, especially with increasing binary eccentricity as also discussed in Sect. 3. Meanwhile, the first stable orbit between the single and the coorbital planet on circular orbits become much more similar with increasing binary eccentricity. In general, we can create a good first order fitting relation by considering a linear function of the first stable pericentre against planet eccentricity as already used in Sect. 3 for the coorbital configuration. The strength of the resonant features will worsen such a fit, especially, when both eccentricities are low for the single planet case.

We based our fit on the quadratic fit function by Holman & Wiegert (1999) and extended it by the planet eccentricity. We reach a best fit using a non-linear least square to the function. For this we take the first stable orbit as a function set and the combination of [e_bin, µ_bin, e_p] as data set to solve for the coefficients assuming a 2nd order term approximation can describe the limits of stability. The mixed terms, e_binµ_bin, e_pµ_bin and $e_{\text{p}}{\mu }_{\text{bin}}^2$ are dropped as they increase complexity without model improvement. Thereby we reach a best fit for the single planet case in the following form: $$\begin{array}{l}{r}_{\text{c},\text{\hspace{0.17em}peri }}/a_{\text{bin}}=(1.36\pm 0.05)\hfill \\ +(5.79\pm 0.19)e_{\text{bin }}-(5.87\pm 0.39)e_{\text{bin }}^2\hfill \\ +(1.99\pm 0.32){\mu }_{\text{bin }}-(3.14\pm 0.52){\mu }_{\text{bin }}^2\hfill \\ +[(1.85\pm 0.05)\hfill \\ -(2.10\pm 0.46)e_{\text{bin }}^2+(3.00\pm 0.52)e_{\text{bin }}{\mu }_{\text{bin }}]e_{\text{p}}.\hfill \end{array}$$(4)

The overall fit reaches a coefficient of determination of R² = 0.962, which is a fair correlation for such a limited data set with three independent variables. However, due to the reduction of the terms the remaining coefficients deviated by factor of at most 4.5 from the values found by Quarles et al. (2018), while describing comparable data. When comparing our zero eccentricity planet case to the study of circular planet orbits by Holman & Wiegert (1999) this leads to a significant deviation, as it ignores the strong resonances acting in circular binaries on circular planets that bend the stability region outwards for planet eccentricity e_p < 0.2 as noticeable in Fig. 8 in the data and a deviation seen in the example fit of Fig. 7. The same resonant lines of instability reaching further pericentre distances have been also observed by e.g. Chavez et al. (2015); Popova & Shevchenko (2016) and Quarles et al. (2018).

In Fig. 9, the first stable coorbital configuration is shown to move further out with increasing binary eccentricity with linear spacing. Changing the binary masses to equal values reduces the difference caused by the binary eccentricity and steepens the slope towards more distant stable orbits for more eccentric planet orbits. However, the stronger non-linear instabilities for e_p > 0.5 make the resulting curve noisy in this region. Nevertheless, from the data we fit a very simple estimate on the stability limits given µ, e_bin and e_p. As our sample size in µ and e_bin is small we used a similarly simple form of equation as Holman & Wiegert (1999) to fit these values to avoid overfitting. We again reduced the equation by all terms that did not produce a contribution above the level of their error after the first fitting iteration. After doing this we arrive at the following simplistic equation to predict an estimate for the first stable orbit for the two-planet coorbital configuration: $$\begin{array}{l}{r}_{\text{c},\text{\hspace{0.17em}peri }}/a_{\text{bin}}=(2.10\pm 0.08)+(2.54\pm 0.12)e_{\text{p}}\hfill \\ +(2.45\pm 0.18)e_{\text{bin}}+(0.90\pm 0.23){\mu }_{\text{bin}}\hfill \\ +(1.37\pm 0.36)e_{\text{p}}{\mu }_{\text{bin}}+(-2.65\pm 0.54)e_{\text{bin }}{\mu }_{\text{bin }}.\hfill \end{array}$$(5)

The fit reaches a coefficient of determination of R² = 0.932 for the limited parameter set in the binary parameters. As the first orbit is nearly linear, this is a good approximation for the realm of stability of eccentric, circumbinary, coorbital configurations. However, as this is just a simple fit, it does not take into account the unstable lines created by the resonances between binary and planet or the unstable regions beyond the first stable orbit for higher planet eccentricities and is only relevant for horseshoe orbit configuration of planets. We investigated some examples of tadpole orbits by initializing the planets near the L4 and L5 points. In the case of tadpole orbits the coorbital planet configuration for all tested parameters grow unstable for e_p > 0.6 as shown in Fig. 10. At this point, all planet configurations with orbit eccentricities e_p > 0.5 should be looked at with caution, while a stability limit for lower eccentricities can be predicted based on our simulations.

Fig. 8 Stability for a parameter range in binary eccentricity (rows) and binary mass ratio (columns) for planet eccentricities against orbit pericentre. Blue indicates unstable orbits, while magenta indicates unstable coorbital conditions. The lines represent the best fit to the first stable orbits.

5 Discussion

Finding horseshoe orbits to be stable is well in agreement with the recent theoretical findings of coorbital multi-planet systems around single stars by Raymond et al. (2023). However, especially the regions of resonance with the binary cause instabilities that extend further for coorbital planet systems than for single planets as the momentum exchange between the planets causes them to enter into the unstable orbit while completing the horseshoe revolution.

The disruption effects of this exchange are comparable to the tidal interaction studied by Mardling (2001). Coorbital resonance is slightly more effective in destabilizing by broadening the effective orbit to the common horseshoe orbit (see e.g. Figs. 3 and 4). This shifts the first stable orbit in the coorbital case compared to the tidal interaction case out by about ~0.1 a_bin for the circular binary. The simulations are limited to two dimensions only. However, most observed, close circumbinary planets are well aligned with the plane of the binary, thereby, our case is the most common picture. The transit observations are biased towards finding primarily well-aligned planets. However, Chen et al. (2020, 2023) showed that the polar planet configuration is more stable than the aligned configuration, therefore our results are conservative cases.

In this study, only the stability of planets that are already initialized in a coorbital configuration is tested without the influence of hydrodynamics. This represents only the final planetary system and does not explain the formation or migration of planets leading to such configurations. The smooth transition from planets embedded in a disc, getting trapped into resonance near the inner cavity, and the disc eventually dissipating would be one of the next steps. In future studies, it would allow us to understand if the 1:1 configuration is affected by the gas dissipation. But the longer simulation time for the hydrodynamic disc model and the n-body is beyond the scope of this work. Martin & Fitzmaurice (2022) investigated the destabilization of single planets during the migration using a simplified static disc model. Their results show paths for stable migration of planets of the mass used in this study to the close-orbits up to the 1:6 resonance with the binary. In Hans Skovgaard), we found using viscous hydrodynamical disc simulations that the inner density maximum close to the inner cavity can provide a high pressure environment which damps migration planets into a 1:1 resonance for nearly circular binaries. _adeg MD) improved the disc description to include eccentricity used in their n-body investigation and found, in fact, some coorbital planets forming. Here, we can show even if the orbits get closer while the disc gets colder and less dense there are still coorbital configurations possible.

However, especially the stable orbits closest in are sensitive to the initial conditions used. Changing the initial conditions such as the precise angle between the pericentre of planets and stars can affect the stability for all simulations of systems with three or more bodies, especially since precession is slow.

For a conservative estimate, the least stable angle was chosen (see also Quarles et al. 2018). While for low eccentricity the first stable orbits still allow for more than 20 full precessions during the simulation time, for e_p = 0.9 the semi-major axis of the first stable orbit is so wide that the planets can only complete about one full precession during the simulation. Thereby the initial configuration can become important. While for a single planet case the long-term stability, especially for these very eccentric planet orbits, could be better determined by a chaos indicator (see e.g. Floud; Sándor et al. 2000), this is not restrictive enough to distinguish between coorbital and non-coorbital two planet systems.

Chavez et al. (2015) and Popova & Shevchenko (2016) investigate the stability for a number of the observed circumbinary planet systems. For low eccentricity binaries (e.g. Kepler-47 in the Popova & Shevchenko 2016 study) or the Kepler-413 case (Chavez et al. 2015; Quarles et al. 2018), the distribution of stable orbits becomes nearly linear especially for e_p < 0.8 and matches well with our finding. Comparing to the classic work of Holman & Wiegert (1999), as we included 3 variable parameters we simplified the fitting function to a lower order. Holman & Wiegert (1999) will be more accurate for circular planet configurations.

When attempting to fit in the same form for the single planet case, we found some deviations with our model, that is strongly informed by eccentric planets. Our fit is generally below their solution by ~0.5 a_bin which also is evident from the top row of Fig. 8. The fit ignores the effects of the resonance. Thus with our e_p-dependent fit, we reach a stability limit that is further in, especially for low planet eccentricities that are disturbed by the most by the resonances. Quarles et al. (2018) gives a good measure for the stability of eccentric single planets taking the strong resonant features into account. Given the focus on eccentric planet orbits in this work, the fit is aimed at a regime where the resonant instabilities lose significance, which is also a natural result for the coorbital planets as they start to become stable further out where the effects of resonances dominate the stability map less.

Previous studies (Quarles et al. 2018; Holman & Wiegert 1999; Popova & Shevchenko 2016; Chen et al. 2019) on single planet stability around binaries have considered a range of binary eccentricities. However, for the 1:1 resonant planet case we put emphasis on nearly circular binaries here, since in Penzlin et al. (2019) we found a hydrodynamical mechanism to form 1:1 resonant planets by migration only around nearly circular binaries. Further, we only consider a limited parameter range to create some general idea about the binary dependent behaviour while rather focussing on the planet eccentricity as a parameter. Fitzmaurice et al. (2022) simulated the migration of multiple planets using rebound with a prescribed disc. Their models achieve coorbital configuration rarely with the circular disc torque used. Such a circular disc torque would in full hydrodynamics simulations occur for more eccentric binaries (Penzlin et al. 2021) than the case studied. However, coorbital resonances are not impossible also for more eccentric host binaries.

Fig. 9 First stable coorbital configurations depending on the planet eccentricity. From left to right the mass ratio µ increases towards equal mass. The colour indicates varying binary eccentricities from circular (blue) to ebin = 0.5 (purple).

Fig. 10 Stability of Tadpole configuration for the Kepler-413 example case. The colour scale indicates the survival time of the planet configuration.

6 Conclusion

In this work, we investigate the stability of coorbital planets compared to single planets around binary stars. In the example of Kepler-413 with one planet, we confirm the previous analysis of circumbinary planet stability (Quarles et al. 2018). In addition, we show that stable coorbital, circumbinary planet configurations are possible even for eccentric planet orbits.

As a good approximation of the 1:1 stability around a Kepler-413 analogue, we found the simple relation of r_c,peri = [(2.90 ± 0.l5)e_p + (2.46 ± 0.08)] a_bin. The stable zone for 1:1 resonant planets is thereby further out than in the single planet case. The critical pericentre for the coorbital configuration is ≥0.8 a_bin farther out than in the case of a single planet, with even greater differences for more eccentric orbits. However, the stability limit of coorbital configuration is still close to the binary and the orbits achieved by migration in Penzlin et al. (2019) are stable for >10⁵ T_bin.

When investigating a range of binary mass ratios and eccentricities, we find that the linear r_c,peri = (A e_p + B) a_bin fit remains a reasonable fit to the data. And with this we can confine stability limits that depend on all three free parameters, µ_bin, e_bin, and e_p, with a very simple first order fit. It is noteworthy that the retrieved stability is specifically valid for coorbital planets in a horseshoe configuration, while tadpole configurations become unstable for planet eccentricity e_p > 0.5. For the horseshoe configuration for e_p > 0.5, unstable regions beyond the first stable orbit appear. For a planet eccentricity e_p < 0.5 and binary eccentricities e_bin < 0.5, we find that close-in, stable, coorbital, planet orbits are possible.

Acknowledgements

In Memory of Willy Kley, who initiated this project. His patient help and kind advice made this possible and we are grateful for every day we worked with him. The authors acknowledge support by the High Performance and Cloud Computing Group at the Zentrum für Datenverarbeitung of the University of Tübingen, the state of Baden-Württemberg through bwHPC and the German Research Foundation (DFG) through grant no INST 37/935-1 FUGG. This research was supported by the Munich Institute for Astro-, Particle and BioPhysics (MIAPbP) which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy -EXC-2094 – 390783311. All plots in this paper were made with the Python library matplotlib (Hunter 2007). A.P. acknowledges support from the Royal Society 2021 Research Fellows Enhanced Research Expenses RF/ERE/210064 and Enhanced Expenses Award and grant KL 650/26 of the German Research Foundation (DFG).

References

All Tables

All Figures

	Fig. 1 Initial configuration of the stars (S1 and S2) and planets (P1 and P2) with the relative argument of periastron Δω between planet and binary orbits.
In the text

	Fig. 2 Stability of coorbital planets with an orbital eccentricity of 0.3 depending on the initial angle between binary periapsis and planet-orbit periapsis. White colour shows stable simulations and the colours indicate the survival time of the coorbital configuration. The dashed and dotted lines mark the 1:x and 2:x planet-binary resonances respectively.
In the text

	Fig. 3 Semi-major axis evolution of the two planets for a time window of 4000 binary orbits after initialization, showing the libration of the planets around a common orbit. The initial semi-major axis was 5.45 abin. However, the average semi-major axis (black dashed line) during the simulation denoted in the top right corner is slightly smaller.
In the text

	Fig. 4 Plot of the relative mean longitude and semi-major axis between the coorbital planets at ap = 5.7 abin with ep = 0.3, indicating that the planets’ relative motion describes a stable horseshoe orbit around the L4 and L5 points marked with red circles, without orbit crossing. The colour progresses over 105 Tbin from yellow to dark blue.
In the text

	Fig. 5 Relative argument of periastron of the common planets orbit and the binary over the precession rate normalized to the theoretical precession rate from the restricted three-body problem.
In the text

	Fig. 6 Map of the survival of single and double (coorbital) planet systems with varying planet eccentricities ep. Blue is unstable for all, purple is only stable for single planets but not the coorbital configuration. The black lines mark l:x resonances between x = 4−13 from shallowest to steepest.
In the text

	Fig. 7 Linear fit of the pericentre of the first stable orbits in Fig. 6. The values close to ep = 0 diverge from the plot due to a unstable region around the 1:4 resonance.
In the text

	Fig. 8 Stability for a parameter range in binary eccentricity (rows) and binary mass ratio (columns) for planet eccentricities against orbit pericentre. Blue indicates unstable orbits, while magenta indicates unstable coorbital conditions. The lines represent the best fit to the first stable orbits.
In the text

	Fig. 9 First stable coorbital configurations depending on the planet eccentricity. From left to right the mass ratio µ increases towards equal mass. The colour indicates varying binary eccentricities from circular (blue) to ebin = 0.5 (purple).
In the text

	Fig. 10 Stability of Tadpole configuration for the Kepler-413 example case. The colour scale indicates the survival time of the planet configuration.
In the text