Issue 
A&A
Volume 680, December 2023



Article Number  A29  
Number of page(s)  8  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202244329  
Published online  05 December 2023 
Stability of coorbital planets around binaries
^{1}
Institut für Astronomie und Astrophysik, Universität Tübingen,
Auf der Morgenstelle 10,
72076
Tübingen, Germany
email: stefan.adelbert@unituebingen.de
^{2}
Astrophysics Group, Department of Physics, Imperial College London,
Prince Consort Rd,
London,
SW7 2AZ,
UK
email: a.penzlin@imperial.ac.uk
^{3}
Department of Applied Mathematics and Physics, Valdosta State University,
Valdosta
GA,
31698, USA
^{4}
Grupo de Dinâmica Orbital e Planetologia, São Paulo State University, UNESP,
Guaratinguetá,
CEP 12516410,
São Paulo, Brazil
Received:
23
June
2022
Accepted:
10
October
2023
In previous hydrodynamical simulations, we found a mechanism for nearly circular binary stars, such as Kepler413, to trap two planets in a stable 1:1 resonance. Therefore, the stability of coorbital configurations becomes a relevant question for planet formation around binary stars. For this work, we investigated the coorbital planet stability using a Kepler413 analogue as an example and then expanded the parameters to study a general nbody stability of planet pairs in eccentric horseshoe orbits around binaries. The stability was tested by evolving the planet orbits for 10^{5} binary periods with varying initial semimajor axes and planet eccentricities. The unstable region of a single circumbinary planet is used as a comparison to the investigated coorbital configurations in this work. We confirm previous findings on the stability of single planets and find a first order linear relation between the orbit eccentricity e_{p} and pericentre to identify stable orbits for various binary configurations. Such a linear relation is also found for the stability of 1:1 resonant planets around binaries. Stable orbits for eccentric horseshoe configurations exist with a pericentre closer than seven binary separations and, in the case of Kepler413, the pericentre of the first stable orbit can be approximated by r_{c,peri} = (2.90 e_{p} + 2.46) a_{bin}.
Key words: binaries: general / planets and satellites: dynamical evolution and stability
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
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^{−5} M_{⋆} for orbit eccentricities up to 0.5. Raymond et al. (2023) recently demonstrated that even multiplanet coorbital configurations can be stable for low masses. Cresswell & Nelson (2006) found in singlestar 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, BalsalobreRuza 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 coorbital 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 planetopened 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^{−5}−10^{−4} 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 nbody simulations. The nbody models, including a prescribed disc by Coleman et al. (2023), also find instances of dynamically created coorbital planet configurations for Kepler16 and Kepler34 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 multiplanets). 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 TOI1338 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 Kepler47 system that hosts three circumbinary planets, demonstrating that multiplanet configurations around binaries are possible (Orosz et al. 2019). In Penzlin et al. (2019), we found stable 1:1 resonant configurations for Kepler47 and 413 analogues in hydrodynamic and nbody 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 zeromass 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 nbody simulations (Rein & Liu 2012; Rein & Spiegel 2015) to identify the stability limit of 1:1resonant, circumbinary planets with a Kepler413like host binary and for a range of binary parameters relevant to the known circumbinary systems.
Section 2 explains the nbody model used to analyse the orbit stability. In Sect. 3, Kepler413 is used as an example to test the best setup, 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 (S_{1} and S_{2}) and planets (P_{1} and P_{2}) with the relative argument of periastron Δω between planet and binary orbits. 
2 Models
In this study, we investigate a binary system hosting circumbinary 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 nbody problem, the variations of parameters, and the simulation setup are defined as follows.
2.1 Binary parameters
We use the Kepler413 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_{2}/M_{bin} = 0.4 using the primary mass M_{1} 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 semimajor 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.
Summary of the physical parameters for the models.
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 coorbital planets are identical, where the first planet starts at the pericentre of the common orbit and the second one at the apocentre. 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 semimajor 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^{−4} 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 planetorbit 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 planetbinary resonances respectively. 
2.3 Simulation setup
To simulate each binary system, we use the rebound nbody code by Rein & Liu (2012) with the IAS15integrator (Rein & Spiegel 2015) for a 2D setup^{1}. To evaluate the stability, we define the following termination conditions: The planetary semimajor 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 semimajor axes is 20% of the inner planet’s semimajor axis evaluated at each time step.
Coorbital configurations must maintain a nearly identical semimajor 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^{5} T_{bin} without termination.
3 Stability of single and 1:1 resonant planets around Kepler413
First, we highlight the stability of single and 1:1 resonant planet systems around a binary similar to Kepler413, because such a loweccentric 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 semimajor axis from the initially set value axis by up to 2% as shown in Fig. 3. The alignment becomes a factor due to the nonperfect 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 binaryplanet resonances, since the libration of the 1:1 resonant planets leads to natural variations in the semimajor 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 semimajor axis as the resonances also become narrow, weakened by the overlap of more higher order resonances and the reduced gravitational potential of the fastmoving 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^{5} T_{bin}. To investigate the dynamics of the planets, we compare our findings with the expected precession period from the restricted threebody 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 threebody 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.5{e}_{\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 threebody problem. This does not change by reducing the planet masses to 10^{−6} M_{bin}.
To compare to previous work we simulate a range of single circumbinary planet systems and vary the initial planetary eccentricity and semimajor axis. Figure 6 shows the instability transition with increasing planetary eccentricity and correspondingly increasing pericentre distance. The nonblue region denotes stable initial conditions for a single planet (up to 10^{5} 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 semimajor 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 Kepler413 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.
Exact Cartesian initial values for all objects in the a_{pini} = 5.45 setup with e_{p} = 0.3 and Δω = π/2.
Fig. 3 Semimajor 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 semimajor axis was 5.45 a_{bin}. However, the average semimajor 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 semimajor axis between the coorbital planets at a_{p} = 5.7 a_{bin} with e_{p} = 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 10^{5} T_{bin} 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 threebody problem. 
Fig. 6 Map of the survival of single and double (coorbital) planet systems with varying planet eccentricities e_{p}. 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 e_{p} = 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 Kepler413like 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 Kepler34 being the most eccentric, planethosting 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^{5} 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 nonlinear 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 \\ \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}+(5.79\pm 0.19){e}_{\text{bin}}(5.87\pm 0.39){e}_{\text{bin}}^{2}\hfill \\ \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}+(1.99\pm 0.32){\mu}_{\text{bin}}(3.14\pm 0.52){\mu}_{\text{bin}}^{2}\hfill \\ \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}+[(1.85\pm 0.05)\hfill \\ \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}(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^{2} = 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 nonlinear 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 twoplanet 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 \\ \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}+(2.45\pm 0.18){e}_{\text{bin}}+(0.90\pm 0.23){\mu}_{\text{bin}}\hfill \\ \text{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}+(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^{2} = 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 multiplanet 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 wellaligned 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 nbody 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 closeorbits 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 nbody 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 semimajor 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 longterm 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 noncoorbital 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. Kepler47 in the Popova & Shevchenko 2016 study) or the Kepler413 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 e_{bin} = 0.5 (purple). 
Fig. 10 Stability of Tadpole configuration for the Kepler413 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 Kepler413 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 Kepler413 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^{5} 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 closein, 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 BadenWürttemberg through bwHPC and the German Research Foundation (DFG) through grant no INST 37/9351 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 EXC2094 – 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
 BalsalobreRuza, O., de GregorioMonsalvo, I., LilloBox, J., et al. 2023, A&A, 675, A172 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chavez, C. E., Georgakarakos, N., Prodan, S., et al. 2015, MNRAS, 446, 1283 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, C., Franchini, A., Lubow, S. H., & Martin, R. G. 2019, MNRAS, 490, 5634 [Google Scholar]
 Chen, C., Lubow, S. H., & Martin, R. G. 2020, MNRAS, 494, 4645 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, C., Lubow, S. H., Martin, R. G., & Nixon, C. J. 2023, MNRAS, 521, 5033 [NASA ADS] [CrossRef] [Google Scholar]
 Coleman, G. A. L., Nelson, R. P., Triaud, A., 2023, MNRAS, 522, 4352 [NASA ADS] [CrossRef] [Google Scholar]
 Cresswell, P., & Nelson, R. P. 2006, A&A, 450, 833 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cresswell, P., & Nelson, R. P. 2009, A&A, 493, 1141 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Doyle, L. R., Carter, J. A., Fabrycky, D. C., et al. 2011, Science, 333, 1602 [NASA ADS] [CrossRef] [Google Scholar]
 Dvorak, R. 1986, A&A, 167, 379 [NASA ADS] [Google Scholar]
 Fitzmaurice, E., Martin, D. V., & Fabrycky, D. C. 2022, MNRAS, 512, 5023 [NASA ADS] [CrossRef] [Google Scholar]
 Froeschlé, C., Gonczi, R., & Lega, E. 1997, Planet. Space Sci., 45, 881 [CrossRef] [Google Scholar]
 Giuppone, C. A., BenítezLlambay, P., & Beaugé, C. 2012, MNRAS, 421, 356 [Google Scholar]
 Holman, M. J., & Wiegert, P. A. 1999, AJ, 117, 621 [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Kley, W., Thun, D., & Penzlin, A. B. T. 2019, A&A, 627, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kostov, V. B., McCullough, P. R., Carter, J. A., et al. 2014, ApJ, 784, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Kostov, V. B., Powell, B. P., Orosz, J. A., et al. 2021, AJ, 162, 234 [NASA ADS] [CrossRef] [Google Scholar]
 Leleu, A., Robutel, P., & Correia, A. C. M. 2018, Celest. Mech. Dyn. Astron., 130, 24 [NASA ADS] [CrossRef] [Google Scholar]
 Leleu, A., Coleman, G., & Ataiee, S. 2019, A&A, 631, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mardling, R. A. 2001, ASP Conf. Ser., 229, 101 [NASA ADS] [Google Scholar]
 Martin, D. V., & Fitzmaurice, E. 2022, MNRAS, 512, 602 [NASA ADS] [CrossRef] [Google Scholar]
 Michalodimitrakis, M., & Grigorelis, F. 1989, J. Astrophys. Astron., 10, 347 [NASA ADS] [CrossRef] [Google Scholar]
 Moriwaki, K., & Nakagawa, Y. 2004, ApJ, 609, 1065 [NASA ADS] [CrossRef] [Google Scholar]
 Muñoz, D. J., & Lithwick, Y. 2020, ApJ, 905, 106 [CrossRef] [Google Scholar]
 Mudryk, L. R., & Wu, Y. 2006, ApJ, 639, 423 [NASA ADS] [CrossRef] [Google Scholar]
 Mutter, M. M., Pierens, A., & Nelson, R. P. 2017, MNRAS, 469, 4504 [NASA ADS] [CrossRef] [Google Scholar]
 Orosz, J. A., Welsh, W. F., Haghighipour, N., et al. 2019, AJ, 157, 174 [NASA ADS] [CrossRef] [Google Scholar]
 Penzlin, A. B. T., Ataiee, S., & Kley, W. 2019, A&A, 630, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Penzlin, A. B. T., Kley, W., & Nelson, R. P. 2021, A&A, 645, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Popova, E. A., & Shevchenko, I. I. 2016, Astron. Lett., 42, 474 [NASA ADS] [CrossRef] [Google Scholar]
 Quarles, B., Satyal, S., Kostov, V., Kaib, N., & Haghighipour, N. 2018, ApJ, 856, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Ragusa, E., Alexander, R., Calcino, J., Hirsh, K., & Price, D. J. 2020, MNRAS, 499, 3362 [NASA ADS] [CrossRef] [Google Scholar]
 Raymond, S. N., Veras, D., Clement, M. S., et al. 2023, MNRAS, 521, 2002 [NASA ADS] [CrossRef] [Google Scholar]
 Rein, H., & Liu, S. F. 2012, A&A, 537, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rein, H., & Spiegel, D. S. 2015, MNRAS, 446, 1424 [Google Scholar]
 Sándor, Z., Érdi, B., & Efthymiopoulos, C. 2000, Celest. Mech. Dyn. Astron., 78, 113 [CrossRef] [Google Scholar]
 Standing, M. R., Sairam, L., Martin, D. V., et al. 2023, Nat. Astron., 7, 702 [NASA ADS] [CrossRef] [Google Scholar]
 Sudol, J. J., & Haghighipour, N. 2021, AJ, 161, 223 [CrossRef] [Google Scholar]
 Sutherland, A. P., & Kratter, K. M. 2019, MNRAS, 487, 3288 [NASA ADS] [CrossRef] [Google Scholar]
 Thun, D., Kley, W., & Picogna, G. 2017, A&A, 604, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Welsh, W. F., Orosz, J. A., Carter, J. A., et al. 2012, Nature, 481, 475 [Google Scholar]
An example setup can be found in github.com/StefanAdelbert/1to1circumbinary
All Tables
Exact Cartesian initial values for all objects in the a_{pini} = 5.45 setup with e_{p} = 0.3 and Δω = π/2.
All Figures
Fig. 1 Initial configuration of the stars (S_{1} and S_{2}) and planets (P_{1} and P_{2}) 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 planetorbit 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 planetbinary resonances respectively. 

In the text 
Fig. 3 Semimajor 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 semimajor axis was 5.45 a_{bin}. However, the average semimajor 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 semimajor axis between the coorbital planets at a_{p} = 5.7 a_{bin} with e_{p} = 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 10^{5} T_{bin} 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 threebody problem. 

In the text 
Fig. 6 Map of the survival of single and double (coorbital) planet systems with varying planet eccentricities e_{p}. 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 e_{p} = 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 e_{bin} = 0.5 (purple). 

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

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.