Primordial migration of co-orbital satellites as a mechanism for the horseshoe orbit of Janus-Epimetheus

TL;DR

This study investigates how primordial migration and disruption of a guiding satellite can naturally lead to the horseshoe co-orbital orbit configuration of Janus and Epimetheus around Saturn, using detailed numerical simulations.

Contribution

It demonstrates that horseshoe co-orbital orbits can emerge naturally from the disruption of a guiding satellite in a circumplanetary disk, providing a new formation mechanism.

Findings

Horseshoe orbits appear as a natural outcome after satellite disruption.

Numerical simulations support the primordial migration hypothesis.

The model explains the current co-orbital configuration of Janus and Epimetheus.

Abstract

We analyze the orbital motion of two natural satellites initially in co-orbital configuration with a third (guiding) satellite embedded into a circumplanetary gas disc and undergoing tidal interactions with the central planet. By solving the exact equations of motion, including the dissipative effects and the mutual gravitational perturbations, we investigate the configuration of the system soon after the guiding satellite is disrupted when crossing the Roche limit with the central planet during its orbital decay. The application is done for the system composed of Saturn and the two small co-orbital satellites Janus and Epimetheus. We perform a large number of numerical simulations, varying the mass of the hypothetical guiding satellite, Saturn's tidal quality factor and the initial configuration and the masses of the small satellites. Analyzing the orbital configurations of Janus and Epimetheus soon after the disruption of the guiding satellite, we obtain that the mutual horseshoe co-orbital orbits appear as a natural outcome of the numerical simulations.

Tables

Table 1. Table 1: Initial conditions used in each set of numerical simulations for the masses of the smallest satellites and their angular configurations. Initially, λ 1 = 0 ∘ subscript 𝜆 1 superscript 0 \lambda_{1}=0^{\circ} for all cases. The mass of the largest satellite and the parameter that controls the orbital decay were changed considering a grid ranging from 0.1 × 10 − 6 ​ m P ≤ m 1 ≤ 1 × 10 − 6 ​ m P 0.1 superscript 10 6 subscript 𝑚 𝑃 subscript 𝑚 1 1 superscript 10 6 subscript 𝑚 𝑃 0.1\times 10^{-6}m_{P}\leq m_{1}\leq 1\times 10^{-6}m_{P} and 0.1 × 10 − 4 ≤ ( k 2 / Q ) P ≤ 1.5 × 10 − 4 0.1 superscript 10 4 subscript subscript 𝑘 2 𝑄 𝑃 1.5 superscript 10 4 0.1\times 10^{-4}\leq(k_{2}/Q)_{P}\leq 1.5\times 10^{-4} . In Simulations 1 to 5, the masses m 2 subscript 𝑚 2 m_{2} and m 3 subscript 𝑚 3 m_{3} correspond to the masses of Janus end Epimetheus, respectively. The last column displays the resulting horseshoe mutual co-orbital configurations after the largest satellite is tidally disrupted as it crosses its Roche limit with the central planet.

Simulation	initial		satellite		horseshoe
	angles		masses		orbits
	${\lambda }_2$	${\lambda }_3$	$m_2$	$m_3$
	(∘)	(∘)	(   ${10}^{-9}{m}_P$)		%
1	60	300	3	1	26.0
2	62	300	3	1	37.8
3	60	298	3	1	19.0
4	62	298	3	1	42.1
5	58	298	3	1	6.73
6	60	300	0.3	0.1	65.0
7	60	300	30	10	0.0

Equations

$${\tau }_{disc}=\frac{a_i}{{\dot{a}}_i}\sim \sqrt{\frac{m_P^3}{\mathcal{G}r}}\frac{h^2}{C_a{m}_i{\sigma }_G},$$

$${\sigma }_G(r)\sim 10g/c{m}^2\sqrt{\frac{1.7R_P}{r}}\frac{0.003}{\alpha }{\frac{0.05}{h}}^2\frac{F_{in}}{10^{-7}{M}_{\odot }y{r}^{-1}}$$

$${\dot{a}}_i=\frac{3n_i{k}_{2P}{m}_i{R}_P^5}{m_P{a}_i^4}{\epsilon }_0,$$

$${\tau }_{tide}=\frac{a_i}{{\dot{a}}_i}=\frac{1}{3n_i}\frac{Q_P}{k_{2P}}\frac{m_P}{m_i}{\left(\frac{a_i}{R_P}\right)}^5,$$

$${\sigma }_G(r)\sim 30g/c{m}^2\sqrt{\frac{1.7R_P}{r}}.$$

Full text

1 Introduction

In the frame of the three-body problem, the co-orbital motion, or 1/1 mean motion resonance (MMR), presents two distinct regimes of libration of the resonant angle, defined by the difference between the mean longitudes of the orbiting bodies, $\Delta\lambda$ (see Murray & Dermott, 1999). When $\Delta\lambda$ librates around 60*∘* or 300*∘* ($L_{4}$ and $L_{5}$ Lagrangian equilibrium solutions, respectively), we have the tadpole regime of motion. When $\Delta\lambda$ librates with high amplitude around 180*∘*, enclosing also $L_{4}$ and $L_{5}$ Lagrangian points, the regime is known as horseshoe motion. In the phase space of the 1/1 MMR, there exist two separatrices for the libration regimes, one separates tadpole and horseshoe orbits and the other one separates the horseshoe orbit from the circulation of $\Delta\lambda$ (see Giuppone et al., 2010 and Rodríguez et al., 2013).

The tadpole orbits are frequent in the Solar System and, in particular, in the saturnian system of natural satellites. Tethys has two small co-orbital satellites, namely, Telesto ($L_{4}$) and Calypso ($L_{5}$), whereas Dione has two others, Helene ($L_{4}$) and Polydeuces ($L_{5}$). Contrarily, the horseshoe orbits are rare: with exception of a few asteroids co-orbital with Earth (De la Fuente Marcos & De la Fuente Marcos, 2016), only Janus and Epimetheus are found in the mutual horseshoe regime of motion (Dermott & Murray, 1981); they are small moons of Saturn orbiting the planet at 2.51 Saturn’s radii on average, sharing the same orbit and swapping their positions every 4 years (Yoder et al., 1983). This difference can be understood from evolutionary models, which predict the formation and dynamical evolution of the systems close to the stationary configurations, particularly, around the points $L_{4}$ and $L_{5}$ (e.g. Jewitt et al., 2004). The horseshoe orbits are, contrarily, peculiar energetic configurations, whose dynamical evolution under external perturbations is easily transformed in a circulation regime of motion (see Rodríguez et al., 2013).

The subject that we are dealing with in this paper is to identify the mechanism that could originate the Janus-Epimetheus mutual horseshoe regime of motion. To our knowledge, up to now, there are no models providing a convincing explanation for the origin of the orbital motion of these small moons, although some works have recently studied the formation of Janus and Epimetheus. Charnoz et al. (2010) and Crida & El Moutamid (2017) suggested that Janus and Epimetheus were formed through aggregation from Saturn’s rings material, whereas Treffenstadt et al. (2015) discussed whether both satellites might have been formed by a collisional disruption of a parental object and evolved into a mutual co-orbital configuration.

Concerning the formation of the saturnian system, Canup & Ward (2006) investigated the possibility that multiple generations of Titan-sized satellites have been lost as a consequence of tidal and gas-driven migration in a primordial circumplanetary disc. Moreover, Canup (2010) argued that the mid-sized saturnian moons, up to Tethys, may have formed by the removal mass from a Titan-sized satellite who crossed its Roche limit. In this context, we purpose the hypothesis that Janus and Epimetheus were co-orbital satellites of an ancient large (guiding) satellite, which was lost during its Roche limit crossing. The fact that Tethys and Dione have their own small co-orbital satellites enforces our hypothesis.

Assuming a primordial circumplanetary disc and tidal interactions, as well as mutual gravitational perturbations, we investigate the orbital configuration of Janus and Epimetheus after the tidal disruption of their hypothetical guiding satellite. We analyze whether both minor satellites might end in a mutual horseshoe regime of motion, testing a large number of initial configurations by varying different parameters involved in the model. For sake of completeness, we also consider more and less massive satellites than Janus and Epimetheus, by an order of magnitude, investigating the implications on the final orbital configuration.

This paper is organized as follows: in Sec. 2 we describe the scenario and the hypothesis adopted in the paper. The interactions considered in the model are detailed in Sec. 3. Sections 4 and 5 display the set of parameters used and the results of the numerical simulations. Finally, Sec.6 is devoted to discussions and conclusions.

2 The scenario adopted

We consider a system composed of a central planet, a main satellite (S1) and two minor satellites (S2 and S3), in such a way that the masses of S2 and S3 are much smaller than the mass of S1. All satellites are initially in co-orbital configuration. We assume that the small satellites are located near the $L_{4}$ (S2) and $L_{5}$ (S3) triangular equilibrium points (see Fig. 1). Following Canup & Ward (2006) and Canup (2010), we adopt a scenario in which the satellites are embedded into a primordial circumplanetary gas disc. This scenario describes the early migration of regular satellites around giant planets corresponding to the final stages of satellite formation. The tidal interaction of the satellite with the central planet is also taken into account, in such a way that the combined effects of disc and tides result in the satellite migration towards the planet (see next section for details).

We assume that the three satellites have an icy composition similar to the mid-sized moons of Saturn (see Table 1 of Canup 2010). The whole system evolves under mutual gravitational perturbation, tidal and disc interactions. We ignore the gravitational influence of companion planets and satellites, restricting the orbital motion to the coplanar case. For the planet-satellite distances that we will consider here, the central star’s gravitational perturbation can also be ignored (see Giuppone and Roig, 2018).

Due to the dissipative forces the satellites undergo inward migration and the Roche limit of S1 ($d_{roche}$) may eventually be crossed, leading to the tidal destruction of the non-differentiated main satellite. We consider that the main satellite is totally disrupted when its Roche limit is attained, while the small satellites continue to evolve under the Keplerian attraction with the central planet and their mutual gravitational perturbation. The aim of this work is to investigate the orbital configuration of the minor satellites when the main migrating satellite crosses its Roche limit.

3 The model

3.1 Gas disc

In the model of Canup & Ward (2006) for gas giant satellite formation, the interaction of the the satellites with the gas disc is modeled by a Type I migration. This type of migration arises due to density waves interactions between the satellite and the circumplanetary gas disc, resulting in orbital decay. In this case, for a satellite at a distance $r$ form the central planet, the timescale of inward migration is given by

[TABLE]

where $a_{i}$ is the semi-major axis of the satellite, $r$ is the instantaneous distance between the planet and the satellite, $m_{P}$ and $m_{i}$ are the masses of the planet and the satellite, respectively, $C_{a}\sim 2.7+1.1q_{G}$, with $q_{G}$ being the surface density profile of the disc ($q_{G}=0.5$, according Canup & Ward 2006), $\sigma_{G}$ is the gas surface density of the circumplanetary disc and $\mathcal{G}$ is the gravitational constant. In addition, $h=c_{s}/\Omega r$ is the vertical height scale or the aspect ratio of the disc, with $\Omega$ being the angular rotation velocity of the local gas and $c_{s}$ is the sound speed. For cold discs, which are rotationally supported, the sound speed $c_{s}$ is smaller than the rotation velocity $v_{rot}=\Omega r$, and $h\sim 0.05$. Following Canup (2010), we assume that Saturn’s rings are in the process of formation, as a consequence of a generation of tidally disrupted moons, and ignore their interaction with the satellites.

In Eq. (1), the Type I migration timescale ($\tau_{disc}$) is inversely proportional to the mass of the satellite ($m_{i}$), hence, smaller satellites migrate slower. In this way, for the masses adopted in this work (see Sec. 4), the rate of orbital decay of the guiding satellite is larger by a factor of 100 when compared to the corresponding rate for the minor satellites. This factor allows us to ignore the influence of the gas disc on the motion of the S2 and S3 satellites.

In our model, $\sigma_{G}$ can be estimated through a quasi steady-state balance between the inflow supply of gas to the disc and the gas removal as the disc viscously spreads into the planet or beyond the outer edge of the disc. In the final stage of the disc evolution, the surface density is given by

[TABLE]

(Canup, 2010), where $R_{p}$ is the planetary radius, $\alpha\sim 0.001$ is the differential rotation of the disc with respect to the Keplerian velocity, $h\sim 0.05$ is the height scale of the disc and $F_{in}$ is the mass inflow rate from the protoplanetary disc, which, for Saturn, is $F_{in}\sim 8\times 10^{-8}M_{\odot}yr^{-1}$ (Canup, 2010).

3.2 Tidal forces

We use the classical Darwin’s formalism to describe the planet-satellite tidal interaction (see Ferraz-Mello et al. 2008, for a review). In this context, the averaged variation of the satellite semi-major axis ($a_{i}$) raised on the planet by the satellite due to tides, is given by

[TABLE]

where $k_{2P}$ is the second-order Love number and radius of the planet and $n_{i}$ is the mean orbital motion of the satellite (above equation holds for small eccentricity, see equation (56) of Ferraz-Mello et al. 2008). $\varepsilon_{0}$ is the tidal phase-lag, which is related to the corresponding tidal frequency $\nu_{o}=2(\omega_{P}-n_{i})$, with $\omega_{P}$ being the angular rotation velocity of the planet. Usually, a linear tidal theory is adopted in such a way that $\varepsilon_{0}=\nu_{0}\Delta t$, where $\Delta t>0$ is the so-called time lag (Mignard 1979; Hut 1981). Hence, the orbital migration due to tides will result in orbital decay or expansion, depending on the sign of $\nu_{0}$. For distances from the planet above the synchronous radius ($\omega_{P}=n_{i}$), we have $\nu_{0}<0$ and, thus, $\dot{a}_{i}<0$, resulting in orbital decay in a timescale given by

[TABLE]

where we have replaced $k_{2P}\,\varepsilon_{0}$ by $k_{2P}/Q_{P}$, with $Q_{P}$ widely used as the tidal quality factor, accounting for the dissipation in the interior of the tidally deformed body.

A simple inspection of Eq. (4) shows that the timescale of orbital decay due to tides in the planet is inversely proportional to $m_{i}$, the mass of satellite, hence, the effect is more important for large satellites. For this reason, we consider only the tidal force acting on the main satellite S1, that corresponds to $i=1$ in Eq. (4).

It is important to note that the semi-major axis will also decrease due to tides raised by the planet on the satellite, however, this contribution is proportional to the square of the eccentricity of the satellite orbit (see Ferraz-Mello et al. 2008) and thus can be neglected for small eccentricities.

4 Numerical simulations

Here we choose Saturn as being the central planet, so $m_{P}=m_{Sat}$, where $m_{Sat}$ is the masss of Saturn. For the minor co-orbital satellites S2 and S3, the masses are $3\times 10^{-9}m_{Sat}$ and $1\times 10^{-9}m_{Sat}$, respectively, which are close to the masses of Janus and Epimetheus111https://ssd.jpl.nasa.gov/horizons.cgi, respectively.

The main satellite S1 is assumed to be homogeneous and have an icy body composition. The Roche limit for an homogeneous (and fluid body) satellite is given by the classical expression $d_{Roche}=2.44R_{P}(\rho_{P}/\rho_{1})^{1/3}$, where $\rho$ stands for mean density (Roche, 1847). As shown in Leinhardt et al. (2012), an homogeneous satellite is fully destroyed when its Roche limit is attained during the inward migration.

On the other hand, in the case of Janus and Epimetheus, the disruption limit is smaller than the classical Roche limit, $d_{Roche}$ (Holsapple and Michel, 2006, 2008; Sharma, 2009). Moreover, according to Porco et al. (2007), there exist the possibility for Janus and Epimetheus to be rubble-pile like objects. Theoretical and numerical works suggest that the Roche limit for a rubble-pile satellite is closer to the planet than for a fluid body of the same mean density (Leinhardt et al., 2012, Zhang et al., 2018). Thus, we assume that S2 and S3 are not disrupted when S1 crosses its Roche limit.

The synchronous radius for the current Saturn’s rotation velocity is $d_{sync}=1.9R_{Sat}$, where $R_{Sat}$ is Saturn’s radius (60,268 km). In addition, replacing values for the current Saturn and an icy satellite, we have $d_{Roche}=2.15R_{Sat}$. Thus, since $d_{sync}<d_{Roche}$, an icy undifferentiated satellite would not undergo inward migration. However, as claimed in Canup (2010), the size and rotation of the early Saturn were different from the current ones, modifying both the synchronous radius and Roche limit values. Following Canup (2010), the synchronous radius was close to 3 – 3.5 $R_{Sat}$ and, adopting $R_{P}=1.5R_{Sat}$ (Canup, 2010), we have $d_{Roche}=2.08R_{Sat}$. Hence, $d_{sync}>d_{Roche}$, allowing for inward migration for $a_{1}\simeq d_{sync}$.

Regarding the expressions for the forces accounting for disc driven and tidal migration, the reader is referred to Tanaka et al. (2002), Tanaka and Ward (2004) and Mignard (1979).

4.1 Initial conditions

For the disc, we choose $\alpha=0.001$ and $h=$0.05, thus the surface density of the gas disc is given by

[TABLE]

Since we consider distances between 2.08 $-$ 2.3 $R_{Sat}$ in our numerical simulations, we have $\sigma_{G}\sim 25g/cm^{2}$.

The mass of S1 is a free parameter in our model and the dissipation within Saturn is not well constrained, thus, we perform several sets of numerical simulations covering a grid in $m_{1}$ and $(k_{2}/Q)_{P}$. The results are shown on a $30\times 20$ grid of initial conditions, allowing us to investigate whether mutual horseshoe co-orbital configurations of S2 and S3 are a possible outcome in the simulations, depending on these parameters.

We vary the mass of S1 (the guiding satellite) in the range (0.1 – 1)$\times 10^{-6}m_{Sat}$. This mass range roughly corresponds to (0.1 – 1)$\times$ Tethys’ mass, one of the icy mid-sized moons of Saturn. We assume a constant density of 1 g/cm3 (around the Tethys’ density), for all considered mass values of S1.

Classical works adopt $Q_{P}>18,000$, based on the assumption that the icy satellites from Mimas to Rhea are primordial in their origin (Goldreich 1965, Peale et al. 1980, Meyer and Wisdom 2007). Recently, evidence for strong tidal dissipation within Saturn was pointed out, equivalent to $(k_{2}/Q)_{P}=1.6\times 10^{-4}$, based on astrometrical measurements of Saturn’s mid-sized icy moons (Lainey et al. 2012, 2017). For Saturn, $k_{2}=0.34$ (Gavrilov and Zharkov, 1977) and the above value of $Q_{P}$ is of the order of 2,000. Therefore, in our numerical simulations, we change $(k_{2}/Q)_{P}$, adopting values in the range (0.1 – 1.5)$\times 10^{-4}$, thus, varying $Q_{P}$ within the range of possible values adopted in classical and recent works.

The initial angular configuration of the system is such that S2 and S3 are placed in the exact triangular equilibrium points, namely, $L_{4}$ and $L_{5}$, with respect to the principal S1 satellite. We set $\lambda_{1}=0$, $\lambda_{2}=60^{\circ}$ and $\lambda_{3}=300^{\circ}$, thus, $\Delta\lambda=\lambda_{3}-\lambda_{2}=240^{\circ}$. We also consider small deviations from $L_{4}$ and $L_{5}$ as initial angular positions of S2 and S3 (see Table 1).

In order to not restrict our results to the specific case of Janus and Epimetheus, we change (increasing and decreasing) the masses of the minor satellites by one order of magnitude.

All simulations start with $a_{1}=a_{2}=a_{3}=2.3R_{Sat}$ (close to the Roche limit for S1, at $2.08R_{Sat}$) and $e_{1}=e_{2}=e_{3}=0.01$, where $e_{i}$ ($i$=1,2,3) are the orbital eccentricities of the satellites. As mentioned above, the early radius of Saturn is $1.5R_{Sat}$ and the initial rotation period is 27.5 h, corresponding to the location of the synchronous radius at $3.5R_{Sat}$ (Canup, 2010). The integration of the equations of motion is performed using a Bulirsch-Stoer integrator with a precision of 10*-13*.

5 Results

The results of the numerical simulations show that, after S1 reaches its Roche limit and is disrupted, the regime of motion of the minor satellites is characterized by the behavior of the angle $\Delta\lambda$. The final results of our simulations can be classified in two configurations:

• •

Regime of mutual tadpole co-orbital motion: $\Delta\lambda$ oscillates around 300∘ (60∘).

• •

Regime of mutual horseshoe co-orbital configuration: $\Delta\lambda$ oscillates with large amplitude around 180∘.

We show in Figs. 2 and 3 examples of possibles final configurations considering two cases with the following initial conditions: $a_{1}$ = 2.3 $R_{Sat}$, $\Delta\lambda=\lambda_{3}-\lambda_{2}=240^{\circ}$, $(k_{2}/Q)_{P}=$ 10*-5* and $\sigma_{G}\sim 25$ g/cm2, only changing the mass of the main satellite: $m_{1}=1\times$ 10${}^{-6}m_{P}$ for the first case (green curves) and $m_{1}=$ 0.9 $\times$ 10${}^{-6}m_{P}$ for the second case (black curves).

Figures 2 and 3 show the temporal evolution of semi-major axes, $\Delta\lambda$ and mean motion ratio between the minor satellites corresponding to the two above examples. The orbital decay leads to the crossing of the Roche limit of S1 ($d_{Roche}$), where the main satellite is disrupted at around 0.325 Myr and 0.35 Myr, for each case. The difference in the timescale for achieving the Roche limit is explained by the difference of the main satellite’s mass ($m_{1}$), that affects the migration velocity (see Eq. (4)). After $d_{Roche}$ is attained, either tadpole (green) or horseshoe (black) co-orbital motions are obtained as final orbital configurations, depending on the mass of S1.

Figure 4 shows, for the same examples, the evolution in the phase space of the 1/1 MMR, where we note the oscillation regimes of the angle $\Delta\lambda$ for both types of co-orbital motions. Red (blue) dots corresponds to the motion before (after) the reaching of $d_{Roche}$. It is interesting to note that, until S1 reaches its Roche limit, $\Delta\lambda$ oscillates with small amplitude around 240*∘*, implying in small oscillations of S2 and S3 around $L_{4}$ and $L_{5}$, respectively, in a tadpole configuration with the guiding body. After the Roche limit, the mutual tadpole (top panel) and horseshoe (bottom panel) co-orbital regimes are obtained in this individual numerical simulations. We note that a small change in the mass $m_{1}$ (from $\times$ 10${}^{-6}m_{P}$ to 0.9 $\times$ 10${}^{-6}m_{P}$) results in two different final orbital configurations of the minor satellites.

The results obtained for the final orbital configurations varying $(k_{2}/Q)_{P}$ and $m_{1}$ according to described in Sec. 4.1, is shown in Fig. 5. We consider $m_{2}=3\times 10^{-9}m_{P}$ and $m_{3}=1\times 10^{-9}m_{P}$ (approximately the masses of Janus and Epimetheus, respectively), whereas $\lambda_{1}=0^{\circ}$, $\lambda_{2}=60^{\circ}$ and $\lambda_{3}=300^{\circ}$, that is, S2 and S3 placed in the exact $L_{4}$ and $L_{5}$ equilibrium points with respect to S1, respectively. This case corresponds to Simulation 1 in Table 1, where the last column indicates the percentage of initial conditions resulting in horseshoe configurations. Green and black squares correspond to tadpoles and horseshoe orbits, respectively. We note that 26% of the initial conditions result in mutual horseshoe orbital configurations between the minor satellites. Moreover, according to Fig. 5, the percentage of final horseshoe orbits decreases for smaller $(k_{2}/Q)_{P}$, indicating that a slow orbital decay scenario (implying in small dissipation within Saturn) is not favorable to explain the occurrence of mutual horseshoe regimes between Janus and Epimetheus.

In order to investigate whether the obtained results are affected by changing the initial angular configuration of the satellites, we perform additional sets of numerical simulations varying the initial mean longitudes, as shown in Table 1 (Simulations 2 to 5). Fig. 6 shows that when S2 is shifted from $L_{4}$, with S3 remaining in $L_{5}$ (Simulation 2, top left panel), the occurrence of horseshoe orbits is almost 1.5 times larger in comparison with the first set of simulations (compare Simulations 1 and 2 in Table 1). On the contrary, changes in the initial position of the minor satellite (S3) does not have significant implications for the final number of horseshoe orbits, as can be seen in Simulations 1 and 3 in Table 1 and the top-right panel in Fig. 6. However, when both minor satellites are moved from $L_{4}$ and $L_{5}$ (Simulations 4 and 5) the results show high and low occurrence of final horseshoe orbits (bottom panels in Fig. 6). In summary, for initial $\lambda_{2}>60^{\circ}$ we obtain the largest number of final horseshoe co-orbital locking in our numerical simulations.

Looking at Figs. 5 and 6 we note that, for high $(k_{2}/Q)_{P}$ (high migration rates), the probability to found horseshoe orbits, as final orbital configuration between the minor satellites, is larger than for small $(k_{2}/Q)_{P}$ values. High value of $(k_{2}/Q)_{P}$ is in agreement with recent determinations for the dissipation within Saturn from astrometric measurements of their natural satellites (e.g., Lainey et al. 2017). On the other hand, the dependence on the mass of the guiding satellite ($m_{1}$) is not clear, although there is a slight trend to get more horseshoe orbits for large $m_{1}$.

We also performed additional sets of numerical simulations modifying the initial positions of the minor satellites, placing the largest one at (or close to, according to the same displacement indicated in Table 1) $L_{5}$ and the smaller one at (or close to) $L_{4}$. The results did not show significant differences with the above sets of numerical simulations shown in Table 1.

As mentioned in Sec. 4.1, we also change the masses of the minor satellites by an order of magnitude smaller and larger than the masses of Janus and Epimetheus (Simulations 6 and 7). We can see form Table 1 and Fig. 7 that the occurrence of horseshoe orbits increases for smaller masses of S2 and S3. Note that the value of $m_{2}/m_{3}$ does not vary in Simulations 1, 6 and 7, indicating that the final orbital configuration in horseshoe-like orbit depends on the individual masses of the minor satellites.

For sake of completeness, we performed additional simulations, now including companion satellites, and analyzed their influence on the final configuration of the system under study. We included the satellites Prometheus and Pandora, which are located close to Janus and Epimetheus, but not Pan, Daphnis and Atlas due to their extremely small masses. We repeated the simulations shown in Fig. 5, showing that there is no significant modifications of the previous percentage of horseshoe and tadpoles final orbits.

6 Discussion and conclusions

In this paper, we analyze the orbital evolution of a system composed of a central planet and three co-orbital satellites. Our model accounts for gravitational interactions between the bodies, a circumplanetary gas disc and tidal interactions. We assume that two (less massive) satellites are initially located close the triangular equilibrium points ($L_{4}$ and $L_{5}$) with respect to the planet and the large (guiding) satellite (Fig. 1). Disc and tidal interactions result in orbital decay of the satellite system, in such a way that the three satellites remain in the co-orbital configuration, with the minor ones in a tadpole regime of motion. Through the solution of the exact equations of motions, we analyze the orbital behavior of the system under assumption that the guiding satellite is completely destroyed soon after it crosses the Roche limit with the planet. Several sets of numerical simulations are performed varying the the mass of the major satellite ($m_{1}$), the migration velocity (parameterized by $(k_{2}/Q)_{P}$) and the initial angular separation of the minor satellites with respect to the exact locations of $L_{4}$ and $L_{5}$ equilibrium points (defined by $\Delta\lambda$).

The results obtained show that, as soon as the guiding satellite crosses its Roche limit, either tadpole or horseshoe regimes of co-orbital motion appear as final configuration between the minor satellites. By varying $m_{1}$ and $(k_{2}/Q)_{P}$ , for a system with Saturn as the central planet and Janus and Epimetheus as the minor satellites, we found that 26% of the numerical simulations result in horseshoe orbits for an initial angular configuration in which Janus and Epimetheus started at the $L_{4}$ and $L_{5}$ equilibrium points. Initial displacements of angular positions from $L_{4}$ or/and $L_{5}$ may significantly change the occurrence of the final horseshoe configurations (see Table 1). Moreover, decreasing or increasing the masses of the minor satellites by one order of magnitude results in drastic changes in the final orbital configurations, with a high occurrence of horseshoe orbits for smallest satellites. The migration rate also plays an important role, since a fast orbital decay favors the occurrence of horseshoe orbits in the numerical simulations.

The results obtained indicate that the currently observed horseshoe configuration between Janus and Epimetheus can be explained in the context of an early migration process interrupted by a cataclysmic event, in which their guiding co-orbital satellite is tidally disrupted via Roche limit crossing. Long-term numerical simulations should be carried out in order to investigate the stability of the obtained horseshoe orbits, however, this is out of the scope of the present paper and can be addressed in a further work.

It is worth noting that, in this work, we do not intend to explain the current orbital configuration of Janus and Epimetheus. Indeed, for sake of simplicity, we do not considered some important effects, such as the interactions with neighbor planets or planetary rings. For instance, Ćuk et al. (2018) have shown that secular resonances with the giant planets may have affected the orbital behaviour of the saturnian satellites in the early Solar System. Moreover, the interaction with Saturn’s rings could induce a transition from the current horseshoe-type lock of Janus and Epimetheus to a tadpole orbit in a timescale of tens of millions years (Lissauer et al., 1985; Caudal, 2013). In addition, it is not yet clear whether Janus and Epimetheus are recent (Charnoz et al. 2010) or primordial (Rossignoli et al., 2019) in their origin. However, although our study does not consider such effects, it points out the possibility to reproduce the mutual co-orbital configuration between km-sized satellites, in the horseshoe regime of motion, through the tidal disruption of a migrating co-orbital large satellite. Thus, we expect that the present work gives an insight into the observed horseshoe orbital confinement of the small saturnian satellites.

Acknowledgements

We acknowledge our referee, Dr. Maryame El Moutamid, for her detailed review and for the helpful suggestions, which allowed us to improve the manuscript. This work was supported by the Brazilian CNPq, CAPES and FAPESP agencies. J. Correa-Otto gratefully acknowledges the financial support by CONICET through PIP 112-201501-00525, and the partial financial support by CICITCA, UNSJ through the projects 21/E1079 (2018-2019) and PROJOVI (2018-2019).