The dynamics of co-orbital giant exomoons – applications for the Kepler-1625 b and Kepler-1708 b satellite systems

ABSTRACT

1 INTRODUCTION

To date, more than 5000 extrasolar planets, exoplanets, have been discovered. This increasing population of bodies presents an outstanding diversity of sizes, masses, and orbital characteristics. Based on the planets of our Solar system, one would expect an abundant population of natural satellites around exoplanets, the exomoons.

Exomoons are yet to be confirmed, but several candidates are reported in the literature (Ben-Jaffel & Ballester 2014; Bennett et al. 2014; Hippke 2015; Lewis et al. 2015; Teachey, Kipping & Schmitt 2018; Heller, Rodenbeck & Bruno 2019; Oza et al. 2019; Fox & Wiegert 2021; Kipping et al. 2022). However, these candidates should be regarded with caution. Some of these satellites are either dynamically unlikely, as is the case for the satellites proposed by Ben-Jaffel & Ballester (2014), which would be orbiting outside the Hill sphere of their host planet, or simply false positives, for example, the six candidates proposed by Fox & Wiegert (2021) that were refuted from both an observational perspective (Kipping 2020) and a stability perspective (Quarles, Li & Rosario-Franco 2020a). In addition, there are satellites proposed based only on indirect effects detected on the host planet (Oza et al. 2019). Among these candidates, the two most promising are Kepler-1625 b-I (Teachey et al. 2018) and Kepler-1708 b-I (Kipping et al. 2022). Both bodies are predicted to be planet-like satellites and were proposed based on transit light curves from the Kepler Space Telescope (Kepler).

The signs of Kepler-1625 b-I were first seen in three transit light curves from Kepler (Teachey et al. 2018) and subsequently identified in data from the Hubble Space Telescope (HST; Teachey & Kipping 2018). However, the only transit from HST did not provide enough information to settle the discussion. Moreover, further analysis of the transits found evidence both in favour (Teachey et al. 2020) and against (Rodenbeck et al. 2018; Heller et al. 2019; Kreidberg, Luger & Bedell 2019) the moon hypothesis, leaving Kepler-1625 b-I as a candidate.

Kepler-1708 b-I was the only emerging exomoon candidate from a survey of 70 transiting cool giant exoplanets compiled by Kipping et al. (2022) using Kepler’s archive. Although there is only one transit exhibiting a faint signal supporting the newest candidate, this evidence checked all the criteria applied by the authors. Recently, Cassese & Kipping (2022) showed that Kepler-1708 b-I is unlikely to be detected by HST, which leaves its status unresolved. Further analysis of the transit presented by Kipping et al. (2022) and more data are needed to validate or refute the candidate.

The lack of confirmed exomoons is due to various reasons, mainly technology limitations. The space telescopes Kepler and CoRot carried expectations on the potential detection of exomoons (Szabó et al. 2006; Simon, Szatmáry & Szabó 2007; Kipping, Fossey & Campanella 2009), but only a few candidates could be proposed using data from these facilities. The HST also could give us hints about exomoons. However, only one transit light curve of HST presented signs of natural satellites (Teachey & Kipping 2018).

In the near future, a new generation of space telescopes is expected to expand the horizon of possibilities for exomoon detection. The recently deployed JWST could provide transit light curves to allow the confirmation of the exomoon candidates (Kipping et al. 2022) in addition to detect new ones (Limbach et al. 2021). The 2026 PLATO mission will search for exoplanets around stars brighter than the ones observed by Kepler (Rauer et al. 2014), thus also increasing the possibility of finding exomoons around close-in giant planets (Heller 2018a; Hippke & Heller 2022).

While waiting for the new telescopes, astronomers have focused on developing techniques to search for signs of exomoons on the available data. Most of the efforts are directed towards the analysis of planetary transits, such that the exomoon’s transit or its indirect effects on the planet’s transit can be identified (Sartoretti & Schneider 1999; Simon et al. 2007; Kipping 2009a, b; Heller 2014; Kipping 2021; Teachey 2021). In addition, Hippke & Heller (2022) released PANDORA, the first publicly available transit-fitting software to search for transits of exomoons, which aims to increase the number of scientists looking for exomoons.

Exomoons can also be addressed from a theoretical point of view. In this way, it is possible to study the stability of exomoons, thus constraining the number, mass, and location of satellites around exoplanets.

Due to observation bias, most of the detected exoplanets are giant planets orbiting close to their host star. Naturally, the first studies about the stability of exomoons were focused on satellites around these planets. Barnes & O’ Brien (2002) applied tidal theory and numerical simulations to constrain the mass of satellites around close-in giant planets that would survive tidal migration and be long-term stable. Their results predicted that Earth-like moons could be stable around Jovian planets depending on the star’s mass and the star–planet separation. However, Barnes & O’ Brien (2002) considered the planet’s inner properties and rotation to be stable for billions of years, which does not correspond to reality. Over long periods of time, planets might shrink (Fortney, Marley & Barnes 2007), while their interior cools and hardens (Guenel, Mathis & Remus 2014). Alvarado-Montes, Zuluaga & Sucerquia (2017), Sucerquia et al. (2019), and Sucerquia et al. (2020) described the tidal migration of exomoons using more robust models. Alvarado-Montes et al. (2017), for example, focused on the radius contraction and internal structure evolution of the planets to constrain the tidal migration of satellites around close-in giant planets. The authors proposed that exomoons exposed to tides would have three fates: (1) fall into the planet; (2) orbital detachment, being ejected from the planet’s Hill sphere; (3) migration to a place of stability, in a quasi-stationary orbit. They found that the inward migration leading to fate (1) is suppressed, such as the exomoons will migrate towards an asymptotic maximum distance leading to fate (3) where the satellites will be stable for long periods of time. Later, Sucerquia et al. (2020) named this distance as satellite tidal orbital parking.

Domingos, Winter & Yokoyama (2006) numerically studied the stability of prograde and retrograde orbits of satellites around giant planets. The authors were able to derive analytical expressions for the critical semimajor axis for the stability of exomoons in both prograde and retrograde motion. The estimations presented by Domingos et al. (2006) were recently updated by Rosario-Franco et al. (2020). From an analytic perspective, Donnison (2010) applied the Hill stability criteria (Murray & Dermott 1999) to place limits on the critical separation between planets and stable satellites for various planet-satellite mass ratios. On the other hand, Namouni (2010) showed that Galilean-like satellites are unlikely to survive the inward migration of their host planet. During an inwards migration, the planet’s Hill radius shrinks, such as the orbits of hypothetical satellites would become more unstable, leading to the ejection of the moons. The results of Namouni (2010) could be related to the lack of satellites detected so far around close-in giant planets.

After the announcement of Kepler-1625 b-I, new models for the stability and formation of exomoons were designed specifically for this candidate. Studies have shown that the origins of this candidate could be explained by both capture (Hamers & Portegies Zwart 2018; Heller 2018b; Hansen 2019) and in situ formation (Moraes & Vieira Neto 2020). However, as these models aimed to reproduce the characteristics of Kepler-1625 b-I as reported by Teachey et al. (2018) and Teachey & Kipping (2018), their applicability to other satellite systems cannot be assured (Kipping et al. 2022).

Given the size, mass, and orbital separation of Kepler-1625 b-I (Teachey et al. 2018; Teachey & Kipping 2018), it is expected that this satellite candidate underwent an intense phase of tidal migration after its formation or capture. Assuming that the host planet is a Jupiter-like body, several studies showed that a Neptune-like satellite would survive the tidal interactions with the host planet and be stable for long periods of time (Rosario-Franco et al. 2020; Tokadjian & Piro 2020; Quarles, Li & Rosario-Franco 2020b). In addition to the satellite candidate’s stability, Kollmeier & Raymond (2019) and Rosario-Franco et al. (2020) showed that theoretical submoons could also be possible around Kepler-1625 b-I. Moreover, Moraes et al. (2022) showed that an extra Earth-like satellite could be stable in the Kepler-1625 b satellite system. The authors explored the regions inside the predicted orbit of Kepler-1625 b-I and found that planetary and satellite tides could stabilize inner satellites, such as these bodies will not fall into the planet or migrate outwards and collide with the satellite candidate. Also, the authors pointed out that the formation of mean motion resonances between the satellites is a dynamical mechanism that secures the stability of both satellites, even when one of the moons has an eccentric orbit.

As it was only recently announced, there are only a few theoretical studies assessing the origins and stability of the exomoon candidate Kepler-1708 b-I. Kipping et al. (2022) studied the post-formation tidal evolution of the candidate, assuming that the satellite was formed in situ at twice the Roche limit. The authors found that once the moon is initially beyond the corotation radius of the system with low eccentricity, the satellite migrates outwards to its proposed position. The authors pointed out that this result does not contribute to determining a possible formation scenario for the candidate, owing to the fact that any model that forms a massive satellite in a tight configuration with the planet could reproduce the tidal outwards migration they presented.

As shown by Moraes et al. (2022), more than one massive moon might be stable in the Kepler-1625 b satellite system. Here, we investigate the possibility of extra satellites being stable in a co-orbital configuration with Kepler-1625 b-I and Kepler-1708 b-I. In our Solar system, several satellites can be found exhibiting co-orbital motion, for example, the Saturnian satellite Tethys and its co-orbital companions Telesto and Calypso, and the classic example of satellites in co-orbital motion, Janus and Epimetheus also in the Saturn system.

Since Kepler-1625 b-I is predicted to be a Neptune-like satellite and Kepler-1708 b-I is expected to be a Super-Earth-like body, we aim to search for stable co-orbital companions that are also planet-like bodies. As shown by Gascheau (1843) and Routh (1874), for the general three-body problem the Lagrangian points L₄ and L₅ can be linearly stable depending on the mass of the bodies.

Before going any further describing our approach to the proposed problem, we must answer the following question: Is it possible to have co-orbital systems composed of planet-like bodies?

As mentioned before, we have co-orbital satellites in our Solar system. However, they are usually formed by two bodies with not comparable masses, with the Janus–Epimetheus system as a remarkable exception, around a planet that is significantly bigger. On the other hand, Kepler-1625 b-I and Kepler-1708 b-I are planet-like satellites, and we are proposing that these bodies have planet-like co-orbital companions. In this way, mechanisms that produce co-orbital planets could be applied to planet-like satellites.

Recently, Long et al. (2022) presented strong evidences that dust could be trapped around the Lagrangian points of a young planet candidate immersed on the LkCa 15 disc using data from the ALMA telescope. This potential discovery was already theoretically predicted in the literature. Montesinos et al. (2020) showed that up to cm-sized particles could be trapped around Lagrangian points in protoplanetary discs. The authors identified that the formation of local vortices at L₄ and L₅ in the early stages of planetary formation are responsible for this dust accumulation. If dust can agglomerate at the Lagrangian points of a giant planet, then the in situ formation of Earth-size co-orbital companions might be possible due to the coagulation of these particles into a single massive rocky object (Chiang & Lithwick 2005; Beaugé et al. 2007; Lyra et al. 2009; Giuppone, Benítez-Llambay & Beaugé 2012). In addition, Laughlin & Chambers (2002) used hydrodynamic simulations to show that a pair of Jupiter-like co-orbital planets could be formed by accretion in a protoplanetary disk, and this 1:1 resonance configuration would be sustainable even after the inwards migration of the planets. Other methods such as pull-down capture of a companion into co-orbital orbit, formation due to direct collision (Chiang & Lithwick 2005), convergent migration of multiple protoplanets (Thommes 2005; Cresswell & Nelson 2006) and gravitational scattering of planetesimals by a protoplanet (Kortenkamp 2005) have been proposed as alternative methods to the in situ formation for the origins of co-orbital planets.

For the formation of co-orbital massive satellites, Moraes & Vieira Neto (2020) showed that in-situ in a massive-solid-enhanced circumplanetary disc could explain the origins of Kepler-1625 b-I. Also, the authors showed that other satellites could form as well. This formation mechanism is compatible with the formation of co-orbital planets if dust could accumulate at the Lagrangian point of the satellite. On the other hand, if Kepler-1625 b-I was captured in a compact configuration, this object could capture pre-existing surviving satellites into co-orbital orbits during its outwards migration phase. Heller (2018b) proposed that Kepler-1625 b-I was part of a binary planetary system before being captured and that the other part of the binary was ejected during the capture phase. It is tempting to suppose that if the binary was captured intact, a co-orbital satellite system could survive. However, the author did not explore this hypothesis, and the successful capture of the complete binary seems unlikely. Another possibility is that the planet-like satellite has its own submoons that ended up being detached from the satellite and later captured into a co-orbital configuration. Even though this is a valid hypothesis, given that submoons are stable only to one-third of the satellite’s Hill radius (Rosario-Franco et al. 2020), here we will be simulating co-orbital companions that are much more massive than the mass limit for submoons proposed by Kollmeier & Raymond (2019) and Rosario-Franco et al. (2020).

In order to study the stability of massive co-orbital satellites in the Kepler-1625 b and Kepler-1708 b systems, we will consider two scenarios. First, the simulated systems will be composed of the planet, the satellite candidate, and the co-orbital companion, neglecting the presence of the star. This assumption considers that Kepler-1625 b-I and Kepler-1708 b-I are both well inside the stability limit for exomoon’s stability (Domingos et al. 2006; Rosario-Franco et al. 2020). In this case, we will be working with a general three-body problem. Subsequently, we will add the star of each system to the simulations. As we will see, co-orbital architectures are sensitive to perturbations, and initial conditions, such as even weaker gravitational effects, could break a once stable configuration, thus justifying the presence of the star. In addition, we want to explore the detectability of systems with co-orbital satellites. In this way, we study the influences of co-orbital satellites over the planet’s transit timing variations (TTVs).

The stability of the co-orbital region of Kepler-1625 b-I and Kepler-1708 b-I will be studied considering co-orbital companions with different: masses, sizes, and initial angular positions. The radius of the bodies is taken into account only to allow collisions. Here, we do not consider any non-gravitational effect.

This paper is organized as follows. In Section 2, we present a review of the physical and orbital characteristics of the bodies forming the Kepler-1625 and Kepler-1708 systems, the initial conditions explored for the extra satellite in each case, and our numerical tools. Then, in Sections 3 and 4, we discuss our results regarding the stability and amplitude of libration of co-orbital satellites, the role of resonances on the surviving of the satellites, how the star of each system influences the stability of the co-orbital region, and an analysis of the planet’s TTVs for the cases where co-orbital satellites are possible. In Section 5, we summarize our results and draw our conclusions.

2 MODEL

In this section, we describe some physical and orbital characteristics of the Kepler-1625 and Kepler-1708 systems, the different systems and initial conditions considered, and present the numerical methods used in this work.

2.1 The Kepler-1625 system

The Kepler-1625 system is composed of a G-type star with approximately 8.7 ± 2.1 Gyr (Teachey & Kipping 2018), which names the system. The star is located in the Cygnus constellation, and it is a solar-mass object (⁠|$M_{\star } \sim 1.079\, {\rm M}_{\odot }$|⁠) with radius |$R_{\star } \sim 1.793\, {\rm R}_{\odot }$| (Mathur et al. 2017).

To date, only one planet has been detected in the system, Kepler-1625b. The planet was first detected via planetary transit in 2015 with data from Kepler (Mullally et al. 2015) and confirmed in 2016 (Morton et al. 2016). Kepler-1625 b has a predicted semimajor axis of a_p ∼ 0.87 au (Morton et al. 2016; Heller 2018b) and its believed to have a coplanar and circular orbit (Teachey et al. 2018). The planet has a radius of R_p = 1.18 Jupiter’s radius (R_J), however, its mass is still not well constrained. Recent photodynamical models showed a distribution of mass peaking at M_p = 3 Jupiter’s masses (M_J) as the most likely value for the mass of Kepler-1625 b (fig. 10 from Teachey et al. 2020).

The exomoon candidate Kepler-1625 b-I was initially predicted to be a Neptune-like body, with a semimajor axis of |$\sim 19.1\, R_p$| and a circular inclined orbit (Teachey et al. 2018). The inclination found for the satellite depends on the detrending method used for transit reduction. Teachey & Kipping (2018) found that the candidate’s inclination will be |$42^{+15\circ }_{-18}$| for linear detrending, |$49^{+21\circ }_{-22}$| for quadratic detrending, and |$43^{+15\circ }_{-19}$| for exponential detrending. In the same work, the author refined the semimajor axis of the satellite, also varying from one data reduction to another, |$45^{+10}_{-5}$|⁠, |$36^{+10}_{-13}$|⁠, and |$42^{+7}_{-4}\, R_p$| for linear, quadratic, and exponential detrending, respectively. Because of these uncertainties, many theoretical studies adopted the canonical value of |$40\, R_p$| for the semimajor axis of the satellite (Hamers & Portegies Zwart 2018; Moraes & Vieira Neto 2020; Moraes et al. 2022; Sucerquia et al. 2022), which agrees with the prediction given by Martin, Fabrycky & Montet (2019). However, other authors used different values for the planet-satellite separation (Tokadjian & Piro 2022).

2.2 The Kepler-1708 system

The Kepler-1708 system is formed by a star, a planet, and the recently proposed exomoon. The star is an F-type, Sun-like object located in the Cygnus constellation. The mass and radius of this body are |$M_{\star } \sim 1.1\, {\rm M}_{\odot }$| and |$R_{\star } \sim 1.1\, {\rm R}_{\odot }$|⁠, respectively (Kipping et al. 2022). The age of the system is estimated to be around 3.16 Gyr.

The planet Kepler-1708b was first detected in 2011 by the Kepler mission, but only recently validated by Kipping et al. (2022). The planet has an estimated semimajor axis of |$1.64^{+0.10}_{-0.10}$| au and is classified as a cool-giant. The planet’s eccentricity is not well determined, and only an upper limit can be found, e_p < 0.4. The same is true for the physical characteristics of the planet. Its mass has an upper bound of |$4.6\, M_J$| and a radius of |$R_p \sim 0.89 \, R_J$| (Kipping et al. 2022). Assuming that the planet has a density similar to Jupiter, Tokadjian & Piro (2022) set the mass of Kepler-1708 b to be |$M_p = 0.81\, M_J$|⁠. The exomoon candidate is predicted to be a Super-Earth-like body with radius |$R_s = 2.61^{+0.42}_{-0.43}\, R_{\oplus }$| and an upper limit for the mass of |$M_s \lt 37\, {\rm M}_{\oplus }$| (Kipping et al. 2022). To find a better estimation for the satellite’s mass, Tokadjian & Piro (2022) set the density of this body to be approximately equal to Neptune’s. In this way, one can find |$M_s = 5\, {\rm M}_{\oplus }$|⁠. Kipping et al. (2022) also found estimations for the orbital radius, |$R_s = 11.7^{+3.8}_{-2.2}$| planetary radius, and inclination, |$9^{+38\circ }_{-45}$|⁠, of the candidate. Despite the uncertainties, if confirmed, the satellite is expected to have low inclination. The motion of Kepler-1708 b-I is predicted to be circular around the host planet.

In Table 1, we present the canonical values for the systems Kepler-1625 and Kepler-1708 adopted in this work.

2.3 Initial conditions

2.3.1 Models without the star – local system

To model our systems, we consider a planetocentric coordinate system, with the planet as the central body. The satellite candidate, Kepler-1625 b-I or Kepler-1708 b-I (hereafter ’primary satellite’), and the extra satellite (hereafter ’co-orbital companion’) are initially in a co-orbital configuration, i.e. the satellites have the same semimajor axis and are only angularly separated.

Moraes et al. (2022) showed that the hypothetical Kepler-1625 b satellite system could be stable with two massive exomoons. The authors drew their conclusions after finding that an extra Earth-like satellite could survive in orbits internal to Kepler-1625 b-I. Here, we will extend this work investigating the possibility of co-orbital satellites in this system and in the Kepler-1708 system. To do so, we will explore a wide range of masses for the co-orbital companion since the stability in the three-body problem involving co-orbital bodies is very sensitive to the mass ratio between the co-orbital pair and the central body.

In our simulations, the mass and radius of the planets and the satellite candidates are taken from Table 1. In all cases, the satellites are in circular and coplanar orbits around the respective planet. The planet–satellite separation is also presented in Table 1.

For the Kepler-1625 system, we consider 18 different types of bodies as co-orbital companions, from Mars-sized to Neptune-sized. We varied the masses of these bodies from M₂ = 0.107 to 17.15 M_⊕, with the intermediate bodies having |$M_2 = i\, {\rm M}_{\oplus }$|⁠, i = 1, ⋅⋅⋅, 16. The radius of the satellites was interpolated using cubic splines taking the values of Mars, Earth, and Neptune as inputs.

Similarly, for the Kepler-1708 system, we have 12 different types of co-orbital companions. The smaller and lighter body is a Mars-like companion, while for the other cases, we consider bodies with M₂ = 0.5 to |$5\, M_{\oplus }$|⁠, with |$\Delta M_2 = 0.5\, M_{\oplus }$| and radius interpolated as before.

Once we set the characteristics of the co-orbital companion, we shall explore the initial angular position related to the primary satellite. Fig. 1 illustrates the initial set-up of our systems. We opted to vary the initial angular separation between the co-orbital satellites from θ_min = 30° to θ_max = 90°, which means we will be exploring the surroundings of the Lagrangian equilibrium point L₄. In preliminary investigations, we found that for angular separations lesser than 30°, the co-orbital structure of the satellites is instantaneously destroyed.

Figure 1.

Illustration of our system’s initial set-up. The planet (black circle) is at the origin of the coordinate system, the primary satellite (blue circle) is placed at a_s from the planet, and the co-orbital companion (open circles) is initially placed at a_s from the planet with an angular separation θ measure from the horizontal line that connects the planet and the primary satellite anticlockwise, from 30° to 90°.

Open in new tabDownload slide

Because of the symmetry of the problem, the results and discussions presented in the following sections are the same for the respective regions around L₅. Thus, we will not analyse the case with the co-orbital companion near L₅.

We consider the satellites to be coplanar to the respective planet and initially in circular orbits. However, Kepler-1625 b-I is thought to be in an inclined configuration (Teachey et al. 2018), but since we are neglecting the presence of the star, we can assume that both satellites are in the same plane as the planet.

2.3.2 Models with the star – complete system

After our initial analysis, we will include the star of the systems and study its gravitation effects on the stability of the co-orbital architectures. For example, Kepler-1625 b has a semimajor axis smaller than 1 au, and it is a three Jupiter’s masses planet. In this case, the gravitational interaction between the planet and the star could generate a non-negligible movement on the centre of mass of the system, which ultimately translates as a wobble on the star. This change in the centre of mass of the system will induce additional movement on the planet and thus jeopardize the fate of co-orbital satellites. Also, by adding the star to our simulations we can access the effects of co-orbital moons on the TTVs of the planets.

The set-ups of the simulations with the star are the same as described before without the star. The planets are considered in circular orbits with the semimajor axis taken from Table 1. Regarding the planet’s inclination, Kepler-1708 b and its satellites are considered coplanar to the star. For the Kepler-1625 system, we study the case with the planet coplanar and with an inclination of 45° relative to the star (other inclinations could be chosen, given the uncertainties on this parameter (Teachey et al. 2018).

2.4 Numerical tools

For our study, we rely on numerical simulations to properly follow the time evolution of these systems since all the bodies involved gravitationally interact with each other.

Our numerical simulations were performed using the IAS15 integration scheme (Rein & Spiegel 2015) implemented in the package POSIDONIUS (Blanco-Cuaresma & Bolmont 2016). POSIDONIUS is often used in problems involving tides or other dissipative effects. However, because of familiarity with this numerical package, we found it convenient to make use of IAS15 written in RUST and opted simply to disable all the dissipative effects implemented in POSIDONIUS, computing only the gravitational interaction between the bodies. For comparison, we randomly chose some of our initial conditions and also simulated with the widely used REBOUND (Rein & Liu 2012) package. The results produced by both packages showed good agreement with compatible computational times.

3 RESULTS FOR THE LOCAL SYSTEM: PLANET-SATELLITE-CO-ORBITAL COMPANION

In this section, we present our results regarding the stability, shape, and amplitude of the co-orbital exomoons’ orbits for the systems Kepler-1625 and Kepler-1708 considering a local system composed of the host planet, the satellite candidate, and the co-orbital companion.

3.1 Stability

First, we present our results for the stability of the systems. We consider a system stable if both satellites are still co-orbitals at the end of the simulations. If the co-orbital configuration is destroyed, we label the system as unstable, regardless of the fate of the satellites after they leave their co-orbital architecture.

As the literature about the dynamics of two massive co-orbital bodies is limited, we will also use results from the co-orbital restricted three-body problem as a first approximation of our findings, such that a comparison between these results can be established.

In Fig. 2, we present our grid of initial conditions (mass of the co-orbital companion versus angular separation) for the system Kepler-1625 (left-hand panel) and Kepler-1708 (right-hand panel), respectively. In green, we have the initial conditions that became stable systems, and in red the unstable conditions. The simulations were carried out for 1 Myr. As one can see, in both cases, the initial conditions for stable systems are around L₄ (θ = 60°), which is an equilibrium point in the restricted three-body problem. The stability of L₄ for the general three-body problem was already predicted, but not shown, by Érdi & Sándor (2005); here, we find that to be true for some cases. Also, there is a pronounced asymmetry for conditions around L₄, especially when the co-orbital companion is less massive. These asymmetries are due to the shape of the co-orbital companion’s orbit, which is more elongated in the opposite direction of the primary satellite. This feature will be explored in detail when we address the shape of the satellites’ orbit.

Figure 2.

Grid with the initial conditions for the Kepler-1625 system (left-hand panel) and Kepler-1708 system (right-hand panel). In green, we have the conditions that end up being stable after 1 Myr, and in red the unstable cases. The horizontal dotted lines mark initial conditions at 60° and 90°.

Open in new tabDownload slide

For all the unstable systems, we found one of the following fates for the satellites: (i) collision between the satellites, (ii) collision between one satellite and the planet, and (iii) ejection of one satellite. Thus, for the local systems, we did not find in our simulations systems where both satellites survived after leaving their co-orbital configuration.

For the Kepler-1625 system, we can see a correlation between the stability and mass of the co-orbital companion (the left-hand panel of Fig. 2). The entire co-orbital region is unstable when the extra satellite has |$M_2 = 5 {\!-\!} 8\, {\rm M}_{\oplus }$| and |$M_2 = 11 {\!-\!} 12\, {\rm M}_{\oplus }$|⁠, and not even the equilibrium point at θ = 60° gave birth to stable systems. For |$M_2 = 9 {\!-\!} 10\, {\rm M}_{\oplus }$|⁠, only a small region close to L₄ is stable. On the other hand, for |$M_2 \ge 13\, {\rm M}_{\oplus }$| the number of stable systems increases, which is counterintuitive. As the mass of the co-orbital companion increases, gravitational interactions between the satellites become stronger. In this way, one would expect the systems’ stability to be compromised. However, we found the opposite. As we increased the mass of the co-orbital companions, more stable systems were found.

The same behaviour is seen for the system Kepler-1708 (the right-hand panel of Fig. 2). The unstable region as a function of the secondary satellite’s mass extends from M₂ = 1.5 to |$2.0\, {\rm M}_{\oplus }$| and from M₂ = 3.0 to |$3.5\, {\rm M}_{\oplus }$|⁠. Same as before, in between these two unstable islands, there is a local stability close to L₄ for the systems with |$M_2 = 2.5\, {\rm M}_{\oplus }$|⁠. Similar to the results for the system Kepler-1625, after the unstable valley, we found more stable conditions as the masses of the secondary satellites are increased.

The above-mentioned results suggest that the instability found for certain values of M₂ is caused by some dynamic effect that depends on the mass of the co-orbital companion. This effect will be explored in the following.

3.2 Resonances of libration frequency

In the restricted three-body problem, the motion of a co-orbital particle about the L₄ of the system is composed by the superposition of two motions (Murray & Dermott 1999). The first motion is a long-period motion of an epicentre librating about the L₄ of the system. Around this epicentre, the co-orbital satellite executes a short-period epicyclic motion, such as the final motion will be the summation of these two movements (figs 3.14 and 3.15 from Murray & Dermott 1999).

Although some similarities between the restricted and general three-body problems are expected, as the co-orbital satellites we are simulating are not particles, they will gravitationally affect the primary satellite, in which both bodies will librate with a particular frequency. To illustrate this behaviour, we show in Fig. 3 the motion of the co-orbital (left-hand panel) and the primary satellite (righ-hand panel) in the frame rotating with the initial circular frequency for 142 yr. For this example, we are considering the Kepler-1625 system, where the co-orbital companion is a Mars-size body and the satellites are initially 48° apart from each other. The same pattern of motion was found for the satellites in the Kepler-1708 system.

Figure 3.

Motion of the co-orbital (left-hand panel) and the primary satellite (right-hand panel) in the Kepler-1625 system. The orbits of the satellites are depicted in the rotating frame for 142 yr. The co-orbital companion has |$0.107\, {\rm M}_{\oplus }$| (Mars-size) and the satellites were initially θ = 48° apart from each other. The colour bar indicates the time. The numbered points show the trajectory sequence of both satellites, such as corresponding positions have the same number.

Open in new tabDownload slide

As one can see, the motions of both satellites depicted in Fig. 3 are similar to the motion of a particle about L₄ in the restricted three-body problem. The epicentre of the co-orbital companion is librating around L₄, while a short-period epicyclic pattern is observed in the loops of the satellite’s trajectory. The same is true for the primary satellite, but in this case, the motion is performed near its initial position. The trajectory of both satellites in the rotating frame is tadpole-like, where the amplitude of the orbits is proportional to the perturbations felt by each satellite.

Similar to the restricted case, here the tadpole-like orbits are more elongated in the opposite direction of the primary satellite. In this way, a greater number of stable systems are expected when we place the co-orbital satellites with θ > 60° than otherwise. This feature explains the asymmetries in the distribution of stable conditions around θ = 60° shown in Fig. 1.

The numbered points in both panels of Fig. 3 represent the position of each satellite at the same time. If we follow these points, we notice that the satellites’ motions are synchronized, such as both satellites crossed their initial semimajor axis at the same time (points of closest and farthest approach). Also, while one satellite is in the inner portion of its orbit (closer to the planet) the other is in the outer portion (farther from the planet) (Fig. 4). At their closest approach, the satellites exchange angular momentum, such as the orbit of the smaller satellite shrinks while the orbit of the bigger satellite expands. The ballet performed by the satellites resembles the motion of Janus and Epimetheus in the Saturn system, but with both satellites having tadpole-like orbits (Epimetheus is in a horseshoe orbit).

Figure 4.

Evolution of the semimajor axis versus time of the satellites for the system Kepler-1625, where the co-orbital companion has |$0.107\, {\rm M}_{\oplus }$| (Mars-size) and the satellites were initially 48° apart from each other. The bottom panel is a zoom on the region with a semimajor axis between 39.95 and |$40.05\, R_p$|⁠. The numbered points correspond to the same points presented in Fig. 3.

Open in new tabDownload slide

Taking the ratio λ_1,2/λ_3,4, if we can find commensurabilities between the frequencies that can be written as the ratio of two integers, we will have resonances of librational frequencies.

One should notice that these frequencies (equations 2 and 3) depend on the mass parameter of the system, such as different resonances will appear only for systems with specific values of μ.

As shown by Érdi & Sándor (2005) (their fig. 6), some librational frequencies, in their study the 2:1 and 3:1, can cause instabilities in co-orbital systems. They showed that around the 2:1 resonance, stable systems are not found for any value of eccentricity of the secondary body, thus creating an instability island for certain values of μ. Also, the authors found that the 3:1 resonance causes the number of stable systems to decrease abruptly. However, they still found stability for cases with lower eccentricities of the secondary body. Other librational frequency resonances can be spotted in their work, for example, the 3:2, but similar to the 3:1 resonance, this commensurability only leads to a decrease in the number of stable systems, implying that stability in this case only can be found if the secondary body does not have an eccentric orbit (e < 0.1).

Some results from the restricted three-body problem are expected to hold in the non-restricted case. In Fig. 2, we have unstable regions for certain masses of the co-orbital companion. These results are similar to the ones found by Érdi & Sándor (2005), thus the nature of the instability could be the same.

To estimate the locations of the librational frequency resonances in our systems, we calculate the ratio of the frequencies of motion, λ_{1, 2}/λ_{3, 4} (equations 2 and 3), using the mass parameter |$\bar{\mu }$| defined in equation (4) for all our systems. We expect that this ratio will give us an approximation for the location of the resonances once the frequencies of motion are valid only for the restricted three-body problem.

In Fig. 5, we show the number of stable systems as a function of the mass parameter |$\bar{\mu }$| for the system Kepler-1625 (purple dotted line) and Kepler-1708 (dark-red dotted line), and the ratio λ_{1, 2}/λ_{3, 4} as a function of |$\bar{\mu }$| (represented by the black curve) which we used to find the libration frequency resonances 2:1 (blue dashed line), 5:3 (red dashed line) and 3:2 (green dashed line). To directly compare our results with Érdi & Sándor (2005), we draw the region in cyan, representing an approximation of the instability region they found.

Figure 5.

Number of stable systems as a function of the mass parameter |$\bar{\mu }$| for the system Kepler-1625 (purple dotted line) and Kepler-1708 (dark-red dotted line). The dots mark the value of |$\bar{\mu }$| as we vary the mass of the co-orbital companion. The black curve is the ratio λ_{1, 2}/λ_{3, 4} as a function of |$\bar{\mu }$|⁠, the vertical black dashed line marks the Gascheau’s criterion limit (⁠|$\bar{\mu } = 0.037$|⁠), the blue, red and green dashed lines denote the location of the 2:1, 5:3 and 3:2 librational frequency resonances for the restricted three-body problem and the region in cyan is an approximation of the region of instability found by Érdi & Sándor (2005) (their fig. 6).

Open in new tabDownload slide

From Fig. 5 one can see that the island of instability found by Érdi & Sándor (2005) was recovered in our simulations. This instability is driven by the 2:1 libration frequency resonance. As we increase the mass of the co-orbital companions, and consequently the value of |$\bar{\mu }$|⁠, we found stable satellites only near L₄. In these cases, the 2:1 resonance is still affecting the co-orbital satellites, but close to L₄, the perturbations are weaker.

Near the 5:3 resonance, we also found only unstable systems. This result agrees with the prediction of Érdi & Sándor (2005), where the authors found only a few stable orbits for satellites in nearly circular motion. We also located the 3:2 libration frequency resonance. For this value of |$\bar{\mu }$|⁠, equivalent to |$M_2 = 16\, M_{\oplus }$| (system Kepler-1625) and |$M_2 \sim 4.5\, M_{\oplus }$| (system Kepler-1708), we found eight and nine initial conditions that returned stable systems for the Kepler-1625 and the Kepler-1708 system, respectively. In all cases, both satellites sustained an almost circular orbit around the planet. These results are in good agreement with the findings of Érdi & Sándor (2005) as well.

3.3 Angular instabilities

In addition to the instabilities that appeared because of the variation in the mass of the co-orbital companion, we detected some smaller islands of instability inside larger islands of stability for certain angular separations. In the Kepler-1625 system, there are two cases of angular instabilities. These two cases appear when we set the co-orbital companion to (i) a Mars-sized body initially with θ = 51° and (ii) a 16-Earth masses body initially with θ = 56°. For case (i), the system becomes unstable after ∼1362 yr when a collision between the satellites takes place. In case (ii), the two satellites also collided, but after only ∼62 yr. For the Kepler-1708 system, when the co-orbital satellites have the same masses (⁠|$5\, {\rm M}_{\oplus }$|⁠) we found instability for two specific angles inside a stability island, θ = 54° and 55°. In both cases, the satellites collided with each other before 2 yr.

Even though the nature of the above-mentioned instabilities is the same, to make the manuscript clearer we will separate the analysis of the Kepler-1625 and Kepler-1708 systems.

3.3.1 Kepler-1625 system

Fig. 6 shows a zoom on the angular separation of the co-orbital satellites around the angles we found peculiar instabilities. For these cases, we performed more simulations using Δθ = 0.1° to identify the extension of the instability islands. As one can see the instabilities are local, with small amplitudes (less than 1°). The structure of these instabilities suggests that initially, we had large stable islands, but due to resonant effects, these islands fragmented into smaller ones. In these cases, the unstable regions are dominated by chaotic motion (see fig. 3 in Liberato & Winter 2020).

Figure 6.

Zoom on the angular separation of the isolated instabilities found in the Kepler-1625 system for a Mars-sized co-orbital companion initially at θ = 51° (top panel) and for a 16-Earth masses co-orbital companion initially at θ = 56° (bottom panel).

Open in new tabDownload slide

One should notice that the librational frequency resonances we studied before cannot be responsible for these local instabilities since the frequencies given by equations (2) and (3) are functions of the mass parameter of the system, while these new features are related to the angular separation of the satellites.

To find the potential resonances acting on θ, we will study the time evolution of the angular separation between the co-orbital satellites. In this way, we can apply a Fast Fourier Transform (FFT) to isolate the dominant frequencies in the time series and investigate if resonances are causing the instabilities we observed.

Fig. 7 shows the magnitude of the FFT associated with the frequencies present in the evolution of the angular separation between the co-orbital satellites.

Figure 7.

Magnitude of the Fast Fourier Transform associated with the frequencies present in the evolution of the angular separation of the satellites in the Kepler-1625 system. Top panel: The co-orbital satellite is a Mars-sized body with an initial angular separation of 51°. Bottom panel: The co-orbital satellite is a 16-Earth masses body with an initial angular separation of 56°.

Open in new tabDownload slide

On the top panel of Fig. 7, we have the FFT analysis of the frequencies of θ for the co-orbital pair formed by the primary satellite and the Mars-sized companion. As one can see, there are two peaks of magnitude around |$2.10\times 10^{-9}\, Hz$| and |$5.30\times 10^{-9}\, Hz$|⁠, representing the two dominant frequencies in θ. Taking the ratio of these two frequencies, we find approximately a 5/2 commensurability, which can be understood as a 5:2 resonance between the libration of the co-orbital satellite about L₄ and the angular motion period of the satellites. This third-order resonance was responsible for increasing the eccentricity of the co-orbital companion. In this way, the orbits of the satellites crossed and led to a collision between the bodies.

The same discussion applies to the co-orbital pair of the primary satellite and the 16-Earth masses companion with initial angular separation θ = 56°. The FFT analysis revealed the frequencies |$3.00\times 10^{-7}\, Hz$| and |$4.50\times 10^{-7}\, Hz$| as the dominant ones on the evolution of θ. Once again, comparing these frequencies, we find a 3:2 commensurability. This particular resonance suddenly increased the eccentricity of both satellites resulting in a collision.

3.3.2 Kepler-1708 system

For the Kepler-1708 system, we found an island of instability inside a greater island of stable initial conditions only for the systems where the satellites have the same masses, |$5\, {\rm M}_{\oplus }$|⁠. In this case, the satellites are stable when their initial angular separation is θ = 53° and for θ = 56°−64°, leaving a gap of unstable conditions for θ = 54° and 55°.

To precisely locate the border of this instability, we follow the approach used for the Kepler-1625 system, first refining the initial conditions between θ = 53° and 55° using Δθ = 0.1° and then applying an FFT analysis for the frequencies of θ. Differently from the previous cases, for the Kepler-1708 system, the amplitude of the angular instability is wider than 1° (Fig. 8) extending from 53.1° to 55.2°. In most cases, the satellites collided with each other within one year, which could point toward the presence of strong resonances acting over the satellites’ angular separation.

Figure 8.

Zoom on the angular separation of the isolated instabilities found in the Kepler-1708 system for a 5-Earth masses co-orbital companion initially located at θ = 54° and 55°.

Open in new tabDownload slide

Fig. 9 shows the FFT analysis of the frequencies in the evolution of θ for the initial separations of 54° (top panel) and 55° (bottom panel), respectively. In both cases, taking the ratio between the dominant frequencies, we find approximately 2, which indicates the proximity of the 2:1 resonance. This first-order resonance abruptly increases the eccentricity of the satellite leading to an almost instantaneous collision.

Figure 9.

Magnitude of the Fast Fourier Transform associated with the frequencies present in the evolution of the angular separation of the satellites in the Kepler-1708 system. Top panel: Initial angular separation of 54°. Bottom panel: Initial angular separation of 56°. In both cases the satellites have the same masses, |$5\, {\rm M}_{\oplus }$|⁠.

Open in new tabDownload slide

3.4 Amplitude of motion

As shown before, the orbits of the primary satellite and its co-orbital companion have a tadpole-like shape. However, we did not address the amplitude of these motions since we only showed one example (a system with a Mars-like co-orbital companion with an initial angular separation of θ = 48° in the Kepler-1625 system).

The amplitude of the satellites’ orbits depends on the magnitude of the perturbations they will receive from the other satellite and the planet. On the other hand, these perturbations are proportional to the mass and initial angular separation of the co-orbital pair. For example, less massive co-orbital companions will be more affected by perturbations from the primary satellite than more massive companions. Also, as the L₄ is an equilibrium point, such as the net force at this point is minimum, at this location, the amplitude of the satellites’ motion is expected to be small.

In Fig. 10, we show the amplitude of motion of the co-orbital companion satellites (Δθ₂) for each initial angular separation for the Kepler-1625 (left-hand panel) and Kepler-1708 (right-hand panel) system, respectively. For better visualization, we only present the cases for θ ≤ 60°, but for greater separations, the behaviour follows the same pattern. The primary satellite always starts at θ = 0°. Thus, its angular displacement is not significant compared to the motion of the co-orbital companion. In this way, we will neglect the analysis of this motion.

Figure 10.

Amplitude of libration of the co-orbital companions (Δθ₂) versus their masses for the systems Kepler-1625 (top panel) and Kepler-1708 (bottom panel). The colour scheme represents the initial θ of each satellite.

Open in new tabDownload slide

As expected, the motions of the co-orbital companions are around L₄ (θ = 60°), although not necessarily symmetrical. Also, as the mass parameter |$\bar{\mu }$| decreases, stable orbits with larger amplitudes of libration are possible (Leleu, Robutel & Correia 2018). In our simulations, we did not find satellites in horseshoe orbits. Roberts (2002) showed throughout analytic calculations that horseshoe configurations in the general three-body problem are stable only for |$\bar{\mu }\le 3\times 10^{-4}$|⁠, which is not the case for our systems. For the restricted three-body problem, numerical simulations found this limit to be μ ≤ 9.5 × 10⁻⁴ (Liberato & Winter 2020).

Stable co-orbital companions initially farther from L₄ presented a larger amplitude of motion because of the perturbations they received from the primary satellite at the moment of their closest approach. If the satellites are too close to each other, this close encounter may result in collisions or ejections of the less massive body. On the other hand, for satellites with initial θ > 60°, the amplitude of libration is similar to the ones presented here.

4 RESULTS FOR THE COMPLETE SYSTEM: STAR-PLANET-SATELLITE-CO-ORBITAL COMPANION

In this section, we study the stability of co-orbital exomoons taking into account the star of the systems. The initial conditions for the planet-satellites systems will be the same as in Section 3, only adding the star as the central body and considering the respective planet with the semimajor axis given in Table 1. The planets will be considered in circular and coplanar orbits. For the system Kepler-1625, the planet is predicted to have a non-neglectable inclination Teachey et al. (2018), in this case, we will also simulate the case with I_p = 45°.

In adding the star, we will study the influence of the satellites on the planet’s motion. In this way, we can find the planet’s TTVs characterized by the presence of co-orbital exomoons.

In the following, we analyse the results for the Kepler-1625 and Kepler-1708 systems, separately.

4.1 Kepler-1625 system

For the Kepler-1625 system considering the star, all the simulated cases ended with unstable systems regardless of the mass of the co-orbital companion, its initial angular position, or the inclination of the planet. In Table 2, we present a summary of the outcomes of the simulations separated by the co-orbital companions’ masses.

Even though the satellites are well inside the Hill radius, |$a_s \sim 0.264\, R_{{\rm H},p}$|⁠, the gravitational perturbations of the star are strong enough to disrupt the initial co-orbital architecture of the satellites. The systems fall apart usually within a few centuries, resulting in collisions between the planet and the satellites, collisions between the satellites, ejections of one of the moons, or the exomoon leaving its planetocentric orbit and assuming an orbit around the star, thus becoming a planet or a ploonet (Sucerquia et al. 2019).

From Table 2, one can see that the most common outcome of our simulations was the collision between the satellites, especially for more massive co-orbital companions. This result is expected in the case of unstable systems since the satellites share the same orbit. Also, we point out the high number of systems ending with the collision between the co-orbital companion and the planet. In this case, the satellite is stripped from its co-orbital orbit, suffers an eccentricity excitation, and collides with the parent body.

To consider the possibility of satellites being ejected from the system, we set the distance of 300 au from the star as the maximum distance a body could reach. Former satellites whose orbits extend beyond this limit are considered ejected from the system. Satellites were ejected only for systems with 5-Earth masses satellites as co-orbital companions. For these systems, the close encounters between satellites are very energetic, resulting in the ejection of one of the satellites, usually the smaller one. Although, we also found cases in which the most massive satellite was the first body ejected from the system.

For the system formed with a 3-Earth masses co-orbital companion initially at θ = 65°, we found that the secondary satellite was ejected from its planetocentric orbit in an inner heliocentric orbit. However, the satellite did not collide with the star and survived as a detached moon, or ploonet as named by Sucerquia et al. (2019). We follow the evolution of this body for 45 thousand years, and the ploonet remained in a stable configuration with the rest of the system. It is out of the scope of this work to investigate in detail the long-term evolution of this body, especially because it represented only |$~0.09{{\ \rm per\ cent}}$| of the total sample. However, as shown in Hansen (2022), unbound moons detached from their planet through tidally driven outward migration are likely to collide with the parental body within millions of years. In this way, the ploonet we found could have the same fate long-term.

4.2 Kepler-1708 system

4.2.1 Stability

Different from the Kepler-1625 system, adding the star of the system Kepler-1708 did not significantly change our results regarding the stability of the co-orbital satellites.

In Table 3, we present a summary of the stable system for the cases with and without the star. As one can see, small differences between the results appear only on the border of the stability regions for some systems, while most of the results from Section 3 are recovered. In this case, we argue that the gravitational influences of the star are minor or even negligible for the Kepler-1708 system.

There are mainly two reasons for the star not being relevant for the stability of co-orbital satellites on the Kepler-1708 system: (i) the star–planet separation; (ii) the position of the satellites related to the planetary Hill radius.

Kepler-1708 b is a cool-giant with a semimajor axis of a_p = 1.64 au. As the gravitational force of a body into another is inversely proportional to the distance between them, the influence of the star over the orbit of the planet, and consequently over its satellites, is less significant than what is expected for a close-in planet. For example, in the Kepler-1625 system, the star–planet distance is 0.87 au and we found that co-orbital satellites are unstable if the star is considered.

The exomoon candidate Kepler-1708 b-I is predicted to have a semimajor axis of |$a_s \sim 0.047\, R_{{\rm H},p}$|⁠. At this distance, the gravitational force of the star is supplanted by the presence of the planet. Thus, the satellites will not feel the presence of the star. For the Kepler-1625 system, the satellites are farther from the planet (⁠|$~0.264\, R_{H,p}$|⁠), while the planet is closer to the star. The combination of these two features jeopardized the presence of massive co-orbital satellites to Kepler-1625 b-I.

4.2.2 TTVs considering co-orbital exomoons

Several authors proposed that TTVs could be used to indirectly detect the presence of more bodies in an exoplanetary system, planets (Holman & Murray 2005; Agol et al. 2005; Nesvorný & Vokrouhlický 2014; Agol & Fabrycky 2018), or exomoons (Sartoretti & Schneider 1999; Szabó et al. 2006; Simon et al. 2007; Kipping 2021). This indirect effect manifests itself as fluctuations in the timing of planetary transits and can be used to infer the presence of planets and moons in systems where at least one planet is known to be transiting.

In this section, we present a TTV analysis for the stable systems with co-orbital exomoons for the Kepler-1708 system.

Our synthetic TTVs will be constructed as follows: (i) we set-up a coplanar system, (x, y), with the star in the origin of the system. The planet is positioned at (a_p, 0), and the satellites are placed around the planet; (ii) we define that the observer is located along the positive direction of the x-axis; (iii) the system will be integrated forward in time, with the planet motion being restricted to the (x, y) plane and anticlockwise; (iv) at each half day interval we will verify if the planet crossed the x-axis from negative y to positive y. If this is the case, we take the time before and after the passage through x and apply a bisection method to precisely find the time of crossing, which we define as a transit time; (v) we stop the simulation when 200 transits are obtained; (vi) at last, we remove the linear trend from our transit times applying a linear least square fit to our data. This process is done for systems considering only Kepler-1708 b-I and systems with the exomoon candidate and its co-orbital companion, such as we can measure the contribution of the co-orbital companion to the planet’s TTV. The simulations are performed using the IAS15 integrator from the REBOUND.¹

In Fig. 11, we compare the TTVs generated by a system with only one moon and a system with a |$5\, {\rm M}_{\oplus }$| co-orbital companion. As one can see, the amplitude of the TTV increased by more than 5 min with the addition of a co-orbital satellite with the same mass as the primary satellite. This variation is significant since the original amplitude of the TTV, considering only one moon, is smaller than 9 min, which is an expected outcome. As shown by Kipping (2009a), the amplitude of the TTV is proportional to the mass of the exomoons. Also, we found that the presence of a co-orbital moon increased the periodicity of the TTVs, allowing us to draw a smoother fit for our data.

Figure 11.

TTVs for a planet with only one moon (purple) and for a planet with a co-orbital pair of satellites (green). The satellites have the same masses, |$5\, {\rm M}_{\oplus }$|⁠, and they are initially 60° apart from each other.

Open in new tabDownload slide

For less massive co-orbital companions, the effects on the planet’s TTV are minor. For example, for a Mars-like companion, the amplitude of the TTV only increased about 0.1 min.

We also investigated the influences of the initial angular separation of the satellites on the planet’s TTVs, but no significant changes were found.

5 CONCLUSIONS

In this work, we studied the stability of co-orbital exomoons using the candidates Kepler-1625 b-I and Kepler-1708 b-I as case studies. The proposed exomoons are predicted to be Super-Earth-like satellites. Thus, we opted to work with massive planet-sized satellites as co-orbital companions. We considered bodies with mass and size varying from Mars-like to a body with the same physical attributes of the respective proposed satellite. The latter being chosen so we will have co-orbital satellites with the same mass and size.

We considered two scenarios, with and without the star as the central body. In the first case, the planet is the central object and the system’s dynamics are modelled by the general three-body problem. Adding the star, we increase the complexity and reality of the problem. The gravitational effects of the star will be of great relevance for the Kepler-1625 system since the host planet, in this case, have a semimajor axis smaller than 1 au. Also, considering the star, we can predict the effects of co-orbital satellites on the planet’s TTV, which is not found in the literature yet.

This work aimed to

• verify the conditions for stability of co-orbital massive exomoons on the Kepler-1625 b and the Kepler-1708 b systems;

• study the role of libration resonances on the stability of the systems;

• understand the correlation between the amplitude of libration of the co-orbital companion satellite and the perturbations from the primary satellite;

• measure the effects of the system’s star over co-orbital exomoons;

• provide planetary TTV profiles for systems with co-orbital exomoons.

In the following, we discuss our results for the systems with and without the star separately.

5.1 Local system: planet–satellite–co-orbital companion

As predicted by Gascheau (1843) and stated later by Érdi & Sándor (2005), the vicinity of L₄ was indeed stable for most of the different co-orbital companions we tested. However, we find that for less massive satellites (Mars-like bodies), this region can extend for initial angular separations between 48° and 84° for the Kepler-1625 system and from 51° to 79° for the Kepler-1708 system. The extension of the stable region slowly decreases as we increased the mass and size of the co-orbital companion until we found stable configurations only near L₄, as initially expected. Also, we found instability islands for specific values of co-orbital satellites’ mass, |$M_2 = 5 {\!-\!} 8\, {\rm M}_{\oplus }$| and |$M_2 = 11 {\!-\! }12\, {\rm M}_{\oplus }$| for the Kepler-1625 system, and |$M_2 = 1.5 {\!-\!} 2.0\, {\rm M}_{\oplus }$| and |$M_2 = 3 {\!-\!} 3.5\, {\rm M}_{\oplus }$| for the Kepler-1708 system. For these systems, not even the initial conditions placed at L₄ survived.

To understand the causes of the instabilities as a function of the mass of the co-orbital companion, we went back to the results of the restricted three-body problem and drew comparisons between our results and the classical ones.

In Érdi & Sándor (2005), the authors mapped the stability of co-orbital systems with different mass parameters under the assumptions of the restricted-three body problem, considering two massive bodies and a particle. They found that librational resonances play a decisive role in raising instabilities depending on the mass parameter of the system. Also, some of these resonances, the 2:1 for example, are so strong that the co-orbital region is unstable for particles, even if the secondary body is in initially circular motion around the primary body.

First of all, we showed that similarly to the restricted case, in our systems, the motion of the co-orbital satellites is described by the superposition of two different motions: a short-period epicyclic motion about the epicentre; and a long-period motion of the epicentre about L₄. Thus, we calculate an approximation for the frequencies of these two motions using equations (2) and (3), derived for the restricted-three body problem, with the mass parameter |$\bar{\mu }$| given by equation (4).

From our analysis, we found that the 2:1 and 5:3 libration resonances may be responsible for the islands of instabilities we detected as we increased the mass of the co-orbital satellites to some specific values. Our results corroborate the findings of Érdi & Sándor (2005) and show that some characteristics of the restricted three-body problem may be valid for the general case. Also, we searched for the location of the 3:2 librational resonance and found that the stable systems at this location have satellites with low eccentricity. The same behaviour was spotted in the restricted case.

In addition to the librational resonances driven by the mass parameter of the system, we found isolated unstable initial conditions located inside islands of stability (⁠|$M_2 = 0.107\, {\rm M}_{\oplus }$| and θ = 51°, and |$M_2 = 15\, {\rm M}_{\oplus }$| and θ = 56° for the Kepler-1625 system, and |$M_2 = 5\, {\rm M}_{\oplus }$| and θ = 54°−55° for the Kepler-1708 system). These instabilities appeared as a function of the initial angular separation of the co-orbital pair.

To find the nature of these two unstable conditions, we generated a time series of the evolution of the angular separation of the satellites and applied an FFT analysis to this series to search for the dominant frequencies of the problem for each case.

For all systems, we found two dominant frequencies. Taking the ratio of the respective frequencies for each system we found some resonances between the motion of the satellites. For the Kepler-1625 system, we found that the instability for the case with |$M_2 = 0.107\, {\rm M}_{\oplus }$| and θ = 51° was generated by a 5:2 resonance between the motion of the satellites, while for the system with |$M_2 = 15\, {\rm M}_{\oplus }$| and θ = 56° we found that a 3:2 resonance was the cause of the instability. In these two situations, the resonances increased the eccentricity of the satellites, which ultimately led to a collision between the co-orbital bodies. The same pattern appeared for the Kepler-1708 system, in this case, a 2:1 resonance was responsible for the instabilities when the co-orbital companion was a |$M_2 = 5\, {\rm M}_{\oplus }$| body initially at θ = 54° and 55°.

Also, we investigated the amplitude of libration of the stable satellites. Here, we gave special attention to the motion of the co-orbital companion, once the movement of the primary satellite in the rotating frame is almost negligible.

The amplitude of libration of the co-orbital satellites is proportional to the mass of the satellite and its initial angular separation. As |$\bar{\mu }$| decreases, stable orbits with larger amplitudes of libration become possible (Leleu et al. 2018). This result was confirmed in our simulations. We found that Mars-like satellites have orbits with larger amplitudes when compared with more massive satellites for the Kepler-1625 system. The surviving co-orbital companions in the Kepler-1708 system presented wider motion amplitude when their masses were |$M_2 = 0.5\, M_{\oplus }$| instead of |$M_2 = 0.107\, M_{\oplus }$|⁠, but given the proximity of these values, our conclusions remain unaffected. Horseshoe orbits are not stable for the values of |$\bar{\mu }$| considered here. Thus, we did not find this configuration in our results.

On the other hand, we found that co-orbital companions initially farther from L₄ are likely to present larger amplitude in their motions. For these cases, the perturbations from the primary satellite and the planet are more pronounced and the amplitude of the tadpole orbit described by the co-orbital companions on the rotating frame increases. If the satellites are angularly closer and suffer a close encounter then the system might become unstable. After reaching instability, the co-orbital satellite may collide with the primary satellite or the planet or even be ejected from the system.

5.2 Complete system: star–planet–satellite–co-orbital companion

The addition of the star proved to be catastrophic for the survival of co-orbital satellites in the Kepler-1625 system. We found that despite the satellites being deep within the Hill radius of the planet, the gravitational influence of the star is still enough to break the co-orbital architectures of the satellites. We found that the most common outcome for the satellites is the collision between them. This is expected given that the satellites initially share the same orbit and experiment close encounters when the initial architecture breaks. In conclusion, massive co-orbital satellites are unlikely on the Kepler-1625 system given the initial conditions we assumed for the planet and satellites.

Our results do not affect previous findings regarding the stability of multiple satellites in the Kepler-1625 system. Here, the initial co-orbital configuration is a major constraint of the problem.

On the other hand, the co-orbital satellites in the Kepler-1708 systems are only marginally disturbed by the addition of the star. This is mainly due to planet–star separation and mass ratio and the initial position of the satellites, inside |$5{{\ \rm per\ cent}}$| of the planet’s Hill radius.

Once we found stable co-orbital satellites for the Kepler-1708 system, we analysed the influences of this type of configuration on the planet’s TTVs. As expected, the amplitude of the TTV increased as the mass of the co-orbital companion increased. For |$M_2 = 5\, {\rm M}_{\oplus }$|⁠, we found an increase of about 5 min on the amplitude of the TTV compared with the case with only one moon. For smaller co-orbital companions, these effects are more subtle. The initial angular position of the co-orbital satellites is not relevant for the TTV produced by the planet.

Finally, the results presented here are dependent on the initial conditions adopted. For the Kepler-1625 system, there are papers showing that the planet–satellite separation can be smaller than the adopted |$40\, R_p$|⁠. This consideration alone will place the satellite in a position where the effects of the star will be less relevant, consequently increasing the possibility of finding stable co-orbital satellites. Moreover, for the local systems, we considered that the planet and the satellites are in the same orbital plane, even though the planet is thought to be inclined regarding its star. When adding the star, we simulated the cases with the planet coplanar and inclined (I_p = 45°), which did not change the scenario of instability for co-orbital satellites. The Kepler-1708 system has even more uncertainties than the previous system. We opted to consider the system’s properties presented in the paper that proposed the exomoon candidate (Kipping et al. 2022) and in the study that explored the tidal evolution of the satellite (Tokadjian & Piro 2022). All in all, more details about these systems are needed to build more accurate models.

ACKNOWLEDGEMENTS

RAM dedicates this paper to his late mentor and friend Willy Kley. The authors thank the anonymous referee for the valuable comments and suggestions that significantly improved this manuscript and Muller Lopes for the help with the TTV analysis. This work was possible thanks to the scholarship granted from the Brazilian Federal Agency for Support and Evaluation of Graduate Education (CAPES), in the scope of the Program CAPES-PrInt, process number 88887.310463/2018-00, Mobility number 88887.583324/2020-00 (RAM). RAM, OCW, and DCM thank the financial support from FAPESP (Grant: 2016/24561-0) and CNPq (Grant: 305210/2018-1). This research was supported by resources supplied by the Center for Scientific Computing (NCC/GridUNESP) of the São Paulo State University (UNESP).

DATA AVAILABILITY

The data underlying this paper will be shared on reasonable request to the corresponding author.

Footnotes

REFERENCES