Stability of coorbital planets around binaries

In previous hydrodynamical simulations, we found a mechanism for nearly circular binary stars, such as Kepler-413, 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 Kepler-413 analogue as an example and then expanded the parameters to study a general n -body 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 semi-major 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 Kepler-413, the pericentre of the first stable orbit can be approximated by r c , peri = (2 . 90 e p + 2 . 46) a bin .


Introduction
Coorbital configurations or 1:1 resonances could thus far 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 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 of 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 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 planet-opened ⋆ passed away 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 allow 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 larger distances 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); Chen et al. (2019) have shown that also eccentric single planets close to the binary can be stable, and 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 , we investigate here 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 section 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 section 4. Finally, we discuss the results in 5 and summarize the findings in 6.

Models
In this study, we investigate a binary system hosting circumbinary planets in an initial coorbital configuration (see Fig. 1).We will 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.

Binary parameters
We will 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 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 semi-major axis 1 a bin for one standard length unit and 1 M bin for the standard mass and vary for the parameter investigation "par".The overall planet mass m p is stated in units of M bin .The planets are initialized with orbital eccentricities e p ranging from 0.0-0.9 with a step width of ∆e p and a varying distance at pericentre with a step width of ∆r peri in units of a bin .The relative argument of periastron between planet and binary orbit ∆ω was set to the least stable value.
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.

Planet parameters
The planetary orbits are initialized relative to the binary center 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 apocenter.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 focused the study on that (see also Sec.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 Section 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 .

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 setup1 .
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 5 T bin without termination.

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 like 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 center 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 center of mass of the entire single planet configuration, but with a true anomaly of π, meaning in its apocenter.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 (2006).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 colors 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.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 a bin .However, the average semi-major axis (black dashed line) during the simulation denoted in the top right corner is slightly smaller.
Varying the initial orbital distance leads to varying orbital stability.For less eccentric planet orbits, the closest distance between planet and binary varies less with ∆ω and, in case of a circular planet orbit, the closest distance remains constant at the difference between planet orbit radius and binary pericenter.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, as the color 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 three-body problem derived from Moriwaki & Nakagawa (2004) T prec = 4 3 (q + 1) 2 q a 3.5 (1 − e 2 p ) 2 1 + 1.5e 2 bin T 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 −6 M bin .
To compare to previous work we simulate a range of single circumbinary planet systems and vary the initial plan-Fig.4: Plot of the relative mean longitude and semi-major 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 color progresses over 10 5 T bin from yellow to dark blue.etary eccentricity and semi-major axis.Fig. 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 5 T bin ).Fig. 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. 7: A 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.
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).Fig. 7 shows the fitted trend of stability for the single planet in blue.Due to resonance some orbits further out are unstable (see Section 4), however, for the more general trend using the first stable orbit is sufficient.
The fitted stability limit is using standard deviation errors.
The simple linear fit describes the stability behaviour of a single planet around a Kepler-413 like binary 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 focusing 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 c,peri = [(2.90± 0.15)e p + (2.46 ± 0.08)] a 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.

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 color 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 Sec. 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 Sec. 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 p µ 2 bin are dropped as they increase complexity without model improvement.Thereby we reach a best fit for the single planet case in the following form: 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 3 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 two-planet coorbital configuration: r c,peri /a bin = (2.10 ± 0.08) + (2.54 ± 0.12)e p (5) + (2.45 ± 0.18)e bin + (0.90 ± 0.23)µ bin + (1.37 ± 0.36)e p µ bin + (−2.65 ± 0.54)e bin µ bin 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.

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.Fig. 3&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. (2020Chen et al. ( , 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 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.  of planets of the mass used in this study to the close-orbits up to the 1:6 resonance with the binary.In Penzlin et al. (2019), 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.Coleman et al. (2023) 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 like 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 allow still 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 longterm stability, especially for these very eccentric planet orbits, could be better determined by a chaos indicator (see e.g.Froeschlé et al. 1997;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., 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 focusing 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.

Conclusion
In this work, we investigate the stability of coorbital planets compared to single planets around binary stars.In the example of a 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.15)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 3 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.
The exact Cartesian initial values for all objects in the a p,ini = 5.45 setup with e p = 0.3 and ∆ω = π/2.

Fig. 2 :
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 color shows stable simulations and the colors indicate the survival time of the coorbital configuration.The dashed and dotted lines mark the 1:x and 2:x planet-binary resonances respectively.

Fig. 5 :
Fig. 5: The 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 :
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 1:x resonances between x=4-13 from shallowest to steepest.

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

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