Moon packing around an Earth-mass planet

ABSTRACT

All four giant planets in the Solar system host systems of multiple moons, whereas the terrestrial planets only host up to two moons. The Earth can capture small asteroids as temporary satellites, which begs the question as to how many moons could stably orbit the Earth, or an Earth-mass exoplanet. We perform a series of N-body simulations of closely spaced equal-mass moons in nested orbits around an Earth-mass planet orbiting a Sun-like star. The innermost moon begins near the host planet’s Roche radius, and the system is packed until the outermost moon begins near the stability limit for single moons. The initial spacing of the moons follows an iterative scheme commonly used for studies of compact planetary systems around single stars. For the three-moon system, we generate MEGNO maps to calculate periodic and chaotic regions and to identify the destabilizing mean motion resonances. Our calculations show that the maximum number of moons depends on the assumed masses of the satellites (Ceres-, Pluto-, and Luna-mass) that could maintain stable orbits in a tightly packed environment. Through our N-body simulations, we find stable configurations for up to 7 ± 1 Ceres-mass, 4 ± 1 Pluto-mass, and 3 ± 1 Luna-mass moons. However, outward tidal migration will likely play a substantial role in the number of moons on stable orbits over the 10 Gyr stellar lifetime of a Sun-like star.

Earth, Moon, planets and satellites: dynamical evolution and stability

1 INTRODUCTION

In our Solar system, most planets contain multiple satellites. Notably, the giant planets in the outer Solar system host multiple-moon systems. The only rocky planets that contain natural satellites are Earth and Mars, with n_moons ≤ 2. Given the discrepancies in total number of moons for the giant and terrestrial planets, it is expected that they experience different formation mechanisms and orbital evolution processes.

Several satellite formation theories have been proposed for both regular and irregular satellites around planets. The commonly accepted mechanism for regular satellite formation around giant planets is accretion from a circumplanetary disc. Canup & Ward (2006) found a planet–disc mass ratio limit of ∼10⁻⁴, which is two to three orders of magnitude smaller than the mass fraction between the Earth and Moon. The results are consistent with simulations of discs with lesser initial mass that have resulted in the formation of satellites that are not massive enough to clear out their orbits, but massive enough to start outward migration due to gravitational interaction with the disc (Hyodo, Ohtsuki & Takeda 2015) and locked in mean motion resonance (MMR).

Interactions between giant planets and circumplanetary discs heavily influence the processes of satellite formation, where the physical parameters that shape disc accretion and evolution have been constrained (Coradini, Magni & Turrini 2010, and references therein). Specifically, the pressure and temperature profiles in the circumplanetary nebulae shaped the chemical gradients in the disc. These chemical gradients set the composition of the satellitesimals, which represent the building blocks of the present regular satellites. Additionally, further studies support the formation of natural satellites around Jupiter and Saturn within the framework of a quasi-steady-state system (Batygin & Morbidelli 2020). A large-scale meridional flow of gas inside the planetary Hill sphere is developed in the later stages of planet formation, which feeds a circumplanetary disc that expels gaseous material back into the parent nebula to maintain equilibrium in the system.

Recent studies have attempted to explain the origin of Galilean moon around Jupiter. Cilibrasi et al. (2021) studied less massive satellite systems through three-dimensional radiative simulations and found that only |$\sim 15{{\ \rm per\ cent}}$| of the resulting population is more massive than the Galilean satellites, which indicates that low rates for tidal migrations and resonant captures are uncommon. Madeira, Izidoro & Giuliatti Winter (2021) reproduced the system of Galilean satellites in a gaseous circumplanetary disc around Jupiter. However, their satellites have moderately eccentric orbits (∼0.1), unlike the current real satellites. They propose a pre-existing resonance between Callisto and Ganymede that was broken over time via divergent migration due to tidal planet–satellite interactions. These same effects further damped the orbital eccentricities of these satellites down to their current values (≲ 0.01).

Satellites can be captured into the gravitational field of a planet resulting in a large semimajor axis, eccentricity, inclination, and in retrograde orbits. These irregular satellites offer an important insight into the formation processes of regular satellites that likely formed in prograde rotating accretion discs. Irregular satellites can be captured by dissociation of a planetesimal binary in the planet’s gravity field (Vokrouhlický, Nesvorný & Levison 2008). For example, Triton was likely captured by this process. Vieira Neto, Winter & Yokoyama (2006) shows that the satellites can be captured in prograde orbits (e.g. Leda, Himalia, Ysithea, and Elara by Jupiter) as a gas giant grows within the planetesimal disc.

Planet-packing studies (Smith & Lissauer 2009; Quarles & Lissauer 2018; Bartram et al. 2021; Lissauer & Gavino 2021) incorporate two or more Earth-mass planets orbiting a Sun-like star in low-eccentricity and low-inclination orbits. Two planets in circular orbits can be Hill stable if the fractional orbital separation is greater than 2.4(μ₁ + μ₂)^1/3, where μ₁ and μ₂ are the planets–Sun mass ratios (Gladman 1993). Also, two small planets are stable if the initial semimajor axis difference (Δ) exceeds 2|$\sqrt{3}$| mutual Hill radii, where systems of more than two planets are stable for |$\Delta \, \gtrsim$| 10 (Chambers, Wetherill & Boss 1996). Similar studies of two-planet systems (equal-mass, coplanar, circular orbits) exhibited stable chaos beyond the |$2\sqrt{3}R_{\rm H}$| separation (Marzari 2014). Closely spaced five-planet systems can have shorter lifetimes when the planetary orbits begin with a non-zero initial eccentricity (e_p = 0.05) in contrast to initially circular orbits (Gratia & Lissauer 2021). Most |$({\sim }72{{\ \rm per\ cent}})$| closely packed five-planet systems with inclined orbits have a prompt collision after their first encounter, where a few |$({\sim }1{{\ \rm per\ cent}})$| survive up to 10^7.5 orbits of the innermost planet (Rice, Rasio & Steffen 2018). Funk et al. (2010) studied hypothetical ultracompact systems of up to 10 planets with Neptune-like masses (17 M_⊕) within 0.26 au and showed that some systems were stable even with a perturbing gas giant between 0.3 and 0.5 au.

The instability in the system arises because the energy and the angular momentum are not conserved due to the perturbations by the additional planet(s). The stability time varies linearly with the initial orbital spacing, and the stability time is significantly higher and can be packed twice as closely together for retrograde planets compared to prograde planets (Smith & Lissauer 2009). On the other hand, planets in circumstellar orbits of binary system require higher spacing than for the planets in single stars (Quarles & Lissauer 2018). To further understand the dynamical stability of multiple moons, we follow the studies done on planet packing in circumstellar orbits in binary star systems due to the similarities in orbital architectures (i.e. natural inner and outer boundaries).

Domingos, Winter & Yokoyama (2006) derived a fitting formula for the stability limit of moons following previous studies of planet stability in binary star systems (Rabl & Dvorak 1988; Holman & Wiegert 1999). The fitting formula defines a critical semimajor axis a_c in units of the planetary Hill radius R_{H, p} as ∼0.5 R_{H, p} or ∼0.9 R_{H, p} for satellites around giant planets in both prograde and retrograde orbits, respectively. Eccentric orbits of either the planet or moon reduce these estimates in a nearly linear manner. Rosario-Franco et al. (2020) clarified the stability limit of a moon in a prograde orbit as a fraction of the Hill radius (0.4|$R_{H_{\rm p}}$|⁠) through a series of N-body simulations that considered a wider range of initial mean anomaly for the satellite. Quarles et al. (2021) revisited the stability limit for retrograde orbits where they showed a limit of |$0.67 R_{\rm H,\, p}$| and identified how the outer stability limit for a putative exomoon in the α Centauri AB system varies due to a forced planetary eccentricity.

Earth has captured small bodies in temporary orbits, where these briefly captured rocks (quasi-satellites) either head into the atmosphere to become meteors (or meteoroids) or orbit the Earth until obtaining the necessary escape velocity to leave Earth’s sphere of influence. The recently detected asteroid CD3, a quasi-satellite that remained in orbit for at least three years, is a prime example of this type of capture. Granvik, Vaubaillon & Jedicke (2012) computed the natural Earth satellite capture probability from the near-Earth object (NEO) population as a function of a NEO’s heliocentric orbital elements. This numerical study included 10 million virtual asteroids and only 18 000 were captured in the Earth orbit. They found that the average captured satellites make approximately three revolutions around the Earth in nine months.

The temporary orbits of quasi-satellites around the Earth and that giant planets naturally contain multiple moons prompt a sensible question as to: ‘How many moons can stably orbit the Earth and how massive can they be?’ In this paper, we perform a series of N-body simulations of closely spaced equal-mass moons in nested orbits around an Earth-mass planet orbiting a Sun-like star to determine the maximum number of moons that could stably orbit the Earth and consider a range of three different prototype masses (Ceres-, Pluto-, and Luna-mass). We use the term Luna to identify a natural satellite that is similar in mass and radius to Earth’s moon.

The methodology of our numerical simulations is presented in Section 2, including the design of the system architectures to simulate the multiple moons (up to nine) in a Sun–Earth system. The results in Section 3 consider Ceres-mass, Luna-mass, and Pluto-mass moons to identify the most stable orbital configuration for the maximum number of moons orbiting the Earth-mass planet. A summary of our results and a discussion of the broader context of multiple-moon system are provided in Sections 3 and 4.

2 METHODOLOGY

Earth and Mars are the only terrestrial planets with moons, where Earth hosts a single moon (the Moon or Luna) and Mars has two moons (Phobos and Deimos). Moon formation is a stochastic process, where the amount of material available largely dictates how many moons could form, but the goal of this work is to find the maximum number of moons that could exist with respect to orbital stability constraints.

2.1 Numerical simulations using rebound

We use the general N-body orbital evolution software rebound (Rein & Liu 2012) to examine the orbital stability of many equal-mass satellites orbiting an Earth-like planet, which in turn orbits a Sun-like star. rebound provides two algorithms (WHFast and IAS15) that are well suited to evolve the hierarchical configuration of stars, planets, and moons. The accuracy of the numerical simulations using each algorithm is not substantially different when the initial time-step is set to 5 per cent of the innermost moon’s initial orbital period and an 11th-order symplectic corrector is used (Rein & Tamayo 2015) to minimize the energy error for WHFast. Therefore, we use the WHFast integrator for the sake of numerical expediency.

Each simulation is evaluated up to 10⁷ times the period of the innermost moon P₁. The time-scale for significant orbital evolution due to tides is much longer and thus we do not consider tidal effects in our N-body simulations. Instead, we use a secular tidal model to evaluate the extent of the moons’ outward migration. An initial configuration is deemed stable, when all the moons initially orbiting the host planet are present at the end of the simulation (⁠|$10^7\, P_1$|⁠). The dynamical time-scale for moon systems is very short, where our time-scale greatly exceeds the secular time-scale for the Sun’s forced eccentricity (<100 yr; Andrade-Ines & Eggl 2017). As a result, systems far from stability boundaries will remain stable for billion-year time-scales. Indeed, our own Moon will evolve on to an unstable orbit (Sasaki, Barnes & O’Brien 2012) eventually, but this time-scale is longer than the main-sequence lifetime of the Sun and renders the issue of stability moot due to the possible engulfment of the Earth–Moon system.

Unstable initial conditions are those that result in a close approach (within the Roche radius) with the host planet, collisions between neighbouring moons or when a moon’s apocentre extends beyond the outer stability limit measured in terms of the planet’s Hill sphere (Q_sat > 0.4 R_{H, p}; Rosario-Franco et al. 2020). Although collisions are possible in our simulations, they are rare and scattering events that transport a moon beyond the outer stability limit represent the vast majority of outcomes. An individual simulation is terminated once an instability occurs, which represents the simulation lifetime t_sim. We scale the simulation lifetime by the orbital period of the innermost moon T₁ to obtain the number of orbits completed by the innermost moon, N₁ = t_sim/T₁. In all of the simulations, the host planet begins on an orbit that is identical to the Sun–Earth system using the JPL Horizons lookup feature of rebound so that a_p ≈ 0.999 au and e_p ≈ 0.0167. As a result, the Sun will perturb each moon’s orbit and lead to a small forced eccentricity (Andrade-Ines & Eggl 2017; Quarles et al. 2021).

2.2 System architecture and formulation of orbital spacing

Although we neglect the long-term orbital effects of tides, tidal forces do place a lower limit on how close a smaller body (e.g. planet or moon) could orbit its parent body (e.g. star or planet). Interior to the Roche limit, the orbiting body gets disintegrated by the tidal force when it overcomes the surface gravity. For a moon with a mass m_sat and radius r_sat orbiting a planet with mass m_p and radius r_p, the Roche radius (via the fluid definition) is given as

$$\begin{eqnarray*} R_{\rm Roche} \approx 2.44 ({m_{\rm p}/m_{\rm sat}})^{1/3}r_{\rm sat}, \end{eqnarray*}$$

(1)

or

$$\begin{eqnarray*} R_{\rm Roche} \approx 2.44 ({\rho _{\rm p}/\rho _{\rm sat}})^{1/3}r_{\rm p}, \end{eqnarray*}$$

(2)

which depends on the bulk density of planet ρ_p and moon ρ_sat through m_p/m_sat = (ρ_sat/ρ_p)(r_p/r_sat)³.

In the three-body problem, the Hill sphere (or radius) defines a region of space where a planet’s gravity dominates the host star’s pull. To first approximation in the planetary eccentricity, the Hill radius for a moon is truncated by the host planet’s pericentre by

$$\begin{eqnarray*} R_{\rm H,p} = a_{\rm p}(1-e_{\rm p})\left(\frac{m_{\rm p}}{3M_\star }\right)^{1/3}, \end{eqnarray*}$$

(3)

which depends on the planet’s semimajor axis a_p, eccentricity e_p, mass m_p, and the host star’s mass M_⋆. In the case of large moons, equation (3) requires modification by replacing the planet mass with |$m_{\rm p}^\prime$|⁠, which is the total mass of i moons added to the planet mass (i.e. |$m_{\rm p}^\prime = m_{\rm p} + {\rm i} m_{\rm sat}$|⁠). The Hill radius is a theoretical point of stability at an instant in time, where many numerical simulations have shown that the outer stability limit actually lies within about half of the Hill radius (Domingos et al. 2006; Rosario-Franco et al. 2020).

Each of the moons begins on a circular, coplanar orbit around an Earth-like planet. The initial semimajor axis of each moon is determined by a unitless spacing parameter β for each simulation. A similar procedure has been used for the study of planet packing around single stars (Chambers et al. 1996; Smith & Lissauer 2009; Obertas, Van Laerhoven & Tamayo 2017) and in binary star systems (Quarles & Lissauer 2018). The spacing parameter β is a normalized separation between the two nearby orbits and uses the mutual Hill radius R_{H, m} between adjacent moons for the normalization. The mutual Hill radius calculated using the total mass that lies interior to ith moon is given by |$\widetilde{M}_i = m_{\rm p} + ({\rm i}-1)m_{{\rm sat},i}$|⁠. The mutual Hill radius for two consecutive moons with mass m_i and m_{i + 1} is defined by

$$\begin{eqnarray*} R_{\rm H,m} = (a_i + a_{i+1})X \,\,\,\, {\rm and} \,\,\,\, X = \frac{1}{2}\left[\frac{m_i+m_{i+1}}{3\widetilde{M}_i}\right]^{1/3}, \end{eqnarray*}$$

(4)

where a recurrence relation defines the semimajor axis of the (i + 1)th outer moon using the semimajor axis of the inner moon a_j, the input parameter β, and X (from equation 4) as

$$\begin{eqnarray*} a_{i+1} = a_i\left(\frac{1+\beta X}{1-\beta X}\right). \end{eqnarray*}$$

(5)

The above equations can be used to describe the mutual Hill radius for a wide range of orbital architectures. In our problem, the m_i + m_{i + 1} factor can be replaced with 2m_sat because we use equal-mass moons.

The initial system setup includes the Sun, Earth, and a single moon with its semimajor axis a₁ at 2 R_Roche. The positions of the subsequent moons (up to nine) are then prescribed using equation (5) for a chosen value of β. Following Smith & Lissauer (2009) and Quarles & Lissauer (2018), the initial mean anomaly is set by using multiples of the golden ratio through 2|$\pi$|iϕ rad = |$360i\, \phi$| degrees, where the golden ratio is |$\phi = (1 + \sqrt{5})/2$|⁠. Using the golden ratio in this context allows us to add the moons to the system so that no pair of moons is initialized at conjunction. This also helps to avoid the MMRs between moons because the first- and second-order MMRs can reduce a system’s lifetime (Quarles & Lissauer 2018).

Fig. 1 illustrates the orbital architecture for eight Ceres-mass moons, including their initial angular position (in Fig. 1a). The colour scheme for the moon index is consistently used throughout this paper. In this colour scheme, the Earth is represented by cyan and the orbital dots (red, blue, green,..., grey) denote each moon in the order from inner to outer. The time variation of a particular moon follows the same colour code in later sections. The inner boundary using Earth’s Roche radius R_Roche and the outer boundary using the stability limit for a single moon (0.4 R_H) are indicated by the dashed orange circles, respectively. The axes’ units in Fig. 1 are converted into Hill radius to provide a logical representation of the Earth’s Hill sphere. Since the Hill radius of any outer moon (compared to an immediate inner moon) is longer, the orbital spacing increases going from the innermost to outermost orbit by a factor of ∼1/3. The displayed orbital spacing and the number of orbits between the stability regions are calculated assuming β = 6.

Figure 1.

Panel (a): schematic of an example orbital architecture showing the initial positions of Ceres-mass moons orbiting an Earth-mass planet. The orbital sizes are determined using the orbital spacing parameter, β, described in Section 2.2, where β = 6. The innermost and outermost dotted circular lines (in orange) denote the inner (Roche radius) and outer stability limit (0.4 R_H). The satellites are colour-coded (red, blue, green,.., grey) from the innermost to outermost moon, where this colour scheme is consistent throughout the paper. The moons’ sizes are not to scale. Panel (b): schematic showing how the spacing parameter β is calculated for the two innermost moons (red and blue). The brown body represents the mutual Hill radius. The Hill radii for m₁, m₂, (m₁ + m₂), and radius of the planet are drawn to scale. The semimajor axes (a₁, a₂, and |$\bar{a}$|⁠) have the same units as the Hill radii, but their scale is broken to minimize the plotting area.

Fig. 1(b) is the projection of the Earth and two innermost moons that are drawn within the grey box in Fig. 1(a). This schematic highlights geometrically how the spacing parameter β distributes two consecutive moons. The Hill radius of each moon individually (R_H1 and R_H2) and their mutual Hill radius (R_{H, m}) are scaled by the planet’s radius R_p. The semimajor axis of each moon has a scale break so that all the bodies can fit on the page. The central body (brown) represents the total moon mass (m₁ + m₂) with the mutual Hill radius R_{H, m}. A similar setup is employed when we consider larger masses for the moons (Pluto-mass and Luna-mass). Note that the mutual Hill radius scales with the assumed mass of the moons, where increasing the moon mass from a Ceres-mass to a Pluto-mass also increases the mutual Hill radius by a factor of (m_Pluto/m_Ceres)^1/3 ≈ 2.4. The mutual Hill radius between Luna-mass moons is ∼4.3 times larger than Ceres-mass moons. As a result, more massive moons will necessarily be limited to smaller values of β so that the outermost moon does not exceed the outer stability limit.

2.3 Initial system parameters

The number of moons that an Earth-like planet can host depends on the assumed satellite masses and their spacing that modulate the gravitational interactions between the satellites. We use three categories (Ceres, Pluto, or Luna) as prototypes for different-sized moons in terms of their mass and radius. The satellite prototypes are used because they represent the most massive object in the asteroid belt (Ceres) and the most massive object in the Kuiper belt (Pluto). The Moon (Luna) is used due to its large relative mass/size among the natural satellites. Phase lag, or constant Q, tidal models suggest that more massive satellites than those we consider can escape an Earth-like planet through outward tidal migration (Quarles et al. 2021, see their fig. 9) and thus we limit our study to planet–satellite mass ratios ≲ 0.02. Probing to smaller masses runs into problems, where we need to consider non-gravitational effects (e.g. Yarkovsky effect or Poynting–Robertson drag). Our study is limited to large and massive objects so that such effects are negligible and can be ignored.

The mean density for all three bodies varies, which results in a difference in the Roche radius and the semimajor axis of the innermost satellite. Table 1 provides the initial values used to define the Roche radius that will scale the innermost satellite orbit. Starting the innermost satellite at the Roche radius could bias our results when that satellite’s eccentricity evolves and its pericentre is raised (i.e. q₁ < R_Roche). We begin the innermost satellite on an initially coplanar, circular orbit with a semimajor axis that is two times the Roche radius (a₁ = 2R_Roche).

Table 1.

Initial satellite parameters (mass, radius, and density) that define the innermost orbit a₁ in terms of R_Roche (see equations 1 and 2) for our N-body simulations. The volumetric mean radius of each satellite type is used, where r_p = 6371 km for the Earth. The period is provided for easier comparisons with other known satellite systems.

Body	Mass	Radius	Density	a1	a1	Period
	(M⊕)	(km)	(g cm−3)	(au)	(R⊕)	(d)
Ceres	0.00015	469.7	2.162	0.000288	6.75498	1.010
Pluto	0.0022	1188	1.854	0.000299	7.01298	1.090
Luna	0.0123	1737.4	3.344	0.000247	5.79333	0.807

To define the first trial value β_min in the orbital spacing of the satellites, we use the dynamical results from previous planet-packing studies (e.g. Gladman 1993; Chambers et al. 1996; Smith & Lissauer 2009; Obertas et al. 2017; Quarles & Lissauer 2018) that show a minimum spacing (⁠|$\beta _{\rm min}=2\sqrt{3}$|⁠), where smaller values are unstable (and chaotic) due to the overlap of first-order MMRs (Wisdom 1980). For β ≳ β_min, there is an expected transition regime up to a critical value β_crit that represents a broad plateau of stable configurations (Obertas et al. 2017; Lissauer & Gavino 2021). The extent of the plateau is expected to change slightly for longer integration time-scale, especially near the MMRs, innermost stable beta, and outermost stable beta due to stochastic encounters. However, this does not exclude the existence of stable conditions within such plateaus. Quarles & Lissauer (2018) showed for binary systems that a maximum value β_max signifies another transition regime, but from stable to unstable due to MMRs with an external perturber. We adapt the results from Quarles & Lissauer (2018) (see their equation 3) to calculate the maximum β through the following:

$$\begin{eqnarray*} \beta _{\rm max} = \left(\frac{(a_N/a_1)^{\frac{1}{N-1}}-1}{(a_N/a_1)^{\frac{1}{N-1}}+1}\right) \left(\frac{12m_{\rm p}}{m_i}\right)^{1/3}, \end{eqnarray*}$$

(6)

which depends on the number of moons N, the innermost orbit a₁, the outermost orbit a_N, the planet’s mass m_p, and each satellite’s mass m_sat. For satellite systems with small N, the value of β_max can be large and thus we only iterate up to β = 10. As the number of moons N increases, we can use equation (6) to verify whether an instability transition occurs. Each of our satellite prototypes (Ceres, Pluto, or Luna) is varied in the initial β starting from |$2\sqrt{3}\approx 3.5$| up to 10 through steps of 0.01.

The spacing parameter β increases the semimajor axis of subsequent moons and the associated orbital periods. As a result, the ratio of orbital periods between a pair of moons can start as a near integer ratio to form an MMR. Although we take steps to avoid MMRs through the initial phasing of the moons, some configurations can still enter into the MMR at least temporarily. To identify the expected locations of MMRs as a function of β (Obertas et al. 2017), we use the semimajor axial ratio a_{i + 1}/a_i from equation (5) and apply Kepler’s third law to get

$$\begin{eqnarray*} \frac{P_{i+1}}{P_i} = \left(\frac{1+\beta X}{1-\beta X}\right)^{3/2}, \end{eqnarray*}$$

(7)

for an adjacent pair of moons. Equation (7) provides the expected location of MMRs assuming that the Sun has a negligible influence on the moons and can be generalized to any pair of moons following the formalism for planet packing (Obertas et al. 2017).

3 RESULTS AND ANALYSIS

3.1 Case study 1: Ceres-mass moons

Using Ceres-mass moons, we perform numerical simulations varying the number of moons (n = 3–9) and their orbital spacing through the spacing parameter (⁠|$2\sqrt{3}\le \beta \le 10$|⁠). Each of the simulations is integrated up to 10⁷ orbits of the innermost moon, which begins at twice the Roche radius. Planet-packing studies using a single star (Smith & Lissauer 2009; Obertas et al. 2017) showed that stability is attained for three to five planets when β ∼ 7–10, where the theoretical minimum stable value β_min is |$2\sqrt{3}$| from Hill stability (Gladman 1993). Table 2 provides estimates for the maximum spacing β_max (see equation 6) for each moon prototype. If we assume a maximum spacing equal to the minimum value for stability in single star systems (β ∼ 7), then we would estimate (using Table 2) a maximum of nine Ceres-mass moons that could be stable. However, we confirm this estimate by performing simulations with 10 moons and find that all the simulations are short-lived. Fig. 1(a) demonstrates how much of the parameter space is filled, where nine Ceres-mass moons with β = 6 nearly reach the outer stability (orange dotted circle).

Table 2.

Values of β_max (see equation 6) when varying the number of moons n using the innermost orbit a₁ and mass m_j for each moon type (see Table 1). The red text marks when |$\beta _{\rm max}\lt 2\sqrt{3}$|⁠.

N	βmax
	mCeres	mPluto	mLuna
3	24.85	10.04	5.96
4	17.76	7.17	4.29
5	13.68	5.51	3.31
6	11.08	4.46	2.69
7	9.30	3.75	2.26
8	8.00	3.22	19.4
9	7.02	2.83	1.71
10	6.25	2.52	1.52

In addition to the simulation lifetime (scaled by T₁), we track each satellite’s maximum eccentricity max e_i. For most of our simulations, each satellite begins on a circular orbit and thus the minimum eccentricity is zero. Fig. 2 shows the maximum eccentricity attained by each moon (Fig. 2a) and the simulation lifetime (Fig. 2b) with respect to the initial orbital spacing parameter β for three Ceres-mass moons. Figs 2(c) and (d) are similar to Figs 2(a) and (b), but evaluate five moons. Similarly, Figs 2(e) and (f) evaluate eight moons. The evaluation of nine moons is not shown here as the system was unstable for less than a million years. The maximum eccentricity is colour-coded (top of the figure) and indicates the index of the moon, where i = 1 and 8 refer to the innermost and outermost moons, respectively.

Figure 2.

Maximum eccentricity (colour-coded) attained for an initial spacing parameter β for three (panel a), five (panel c), and eight (panel e) Ceres-mass moons. The system lifetime measured in terms of the number of orbits N₁ completed by the innermost moon is plotted on a logarithmic scale for three (panel b), five (panel d), and eight (panel f) Ceres-mass moons. The orange dashed line marks the maximum β for the highest number of moons (see Table 2). The cyan line shows the log-linear fit to the unstable lifetimes, which is discussed in Section 3.4.

Since the Earth-mass host planet begins with a non-zero eccentricity (e_p ≈ 0.0167), then each satellite experiences a forced eccentricity, which arises in solving the secular quadrupole problem of the orbit-averaged disturbing function (Andrade-Ines & Eggl 2017). The magnitude of the forced eccentricity due to the Sun is typically small, but the eccentricity growth from moon–moon interactions can quickly drive up a moon’s eccentricity enough (e_i ∼ 0.1) for orbit crossings to occur. In Fig. 2(a), the maximum eccentricity of each moon e_i is high due to the chaotic overlap of first-order MMRs (Deck, Payne & Holman 2013) for |$2\sqrt{3}\le \beta \lesssim 6$|⁠. Once β > 6, the gravitational perturbations from moon–moon interactions weaken and the moon pairs exit the chaotic zone. Beyond β ∼ 6, the maximum eccentricity of each moon remains low (⁠|${\lesssim 0.1}$|⁠), but non-zero. There are spikes in the maximum eccentricity (in Fig. 2a) that correlate with dips in the lifetimes |$(\log _{10}\, N_1)$| measured in Fig. 2(b). These anomalous features (⁠|$\beta = 7.2,\, 7.4,\, 7.6,\, \ldots$|⁠) correspond to first-order MMRs between adjacent pairs of moons. The rise in maximum eccentricity at β ∼ 9.8 corresponds to the 2:1 MMR, which is expected at β ≈ 9.78 (Obertas et al. 2017). The stable plateau for three Ceres-mass moons continues until β ∼ 24, where the apocentre of the outermost moon will extend beyond the outer stability limit.

Considering systems with additional satellites modifies Figs 2(a) and (b) by lowering β_max (see Table 2), and the limit on the number of moons is restricted by the minimum value, |$2\sqrt{3}$|⁠. The maximum eccentricities for systems of five and eight Ceres-mass moons are given in Figs 2(c) and (e), respectively, while the lifetimes are shown in Figs 2(d) and (f), respectively. These panels illustrate similar features as Figs 2(a) and (b), where the tail of stability is apparent at β ≈ 8 (see the vertical dashed line in Fig. 2f). For 6.6 ≲ β ≲ 7.6, a system of eight moons is on the border of stability as evidenced by the growth of eccentricity with β for the outermost moon and additional effects (e.g. tides) may prevent long-term stability in general. A dip appears at β ∼ 7.3 in Fig. 2(f), which corresponds to the 5:3 MMR between adjacent moons. The width of the MMRs grows as moons are added to the system because there are more small perturbations possible that can push a moon into an MMR as the system evolves (i.e. evolution of a moon’s semimajor axis or β).

Tracking the maximum eccentricity of each moon’s orbit shows which values of β are likely to produce unstable systems. The time-scale for 10⁷ orbits of the innermost moon is only ∼20 000 yr, where longer simulations are impractical due to the small time-step required for accurately evolving each system. The gravitational interactions between moons also occur on a short time-scale, which further necessitates a small integration time-step. Despite these limitations, the maximum eccentricity provides a good heuristic to show the likely long-term stability of moons. Although not shown here, we perform simulations with four, six, and seven moons and confirm the apparent trends in stability.

3.2 Case study 2: Pluto-mass moons

The mutual Hill radius increases with Pluto-mass moons, which corresponds to a larger physical spacing between moons through the parameter β. The expected locations of the MMRs between moons more strongly depend on β, where they move to lower values by a factor of ∼2.44, or (m_Pluto/m_Ceres)^1/3. With these expectations, we perform simulations of Pluto-mass moons following the same procedure from Section 3.1.

Fig. 3 illustrates how the maximum eccentricity of each moon and system lifetime varies with respect to the spacing parameter β for a system of three (Figs 3a and b) and five (Figs 3c and d) moons. From Table 2, we expect a system of three moons to be stable up to β ∼ 10, where Fig. 3(a) shows the increasing maximum eccentricity attained by the outermost moon (i = 3) as β increases. The minimum spacing appears at β ∼ 4.5 because the 2:1 MMR is expected at β ≈ 4 and planet-packing studies have shown broad stability occurring beyond this MMR (Smith & Lissauer 2009; Obertas et al. 2017; Quarles & Lissauer 2018; Lissauer & Gavino 2021). The underlying physical cause is the wider separation of libration zones for the first-order MMRs and the wider physical separation (i.e. less gravitational perturbations) between moons. The spike in Fig. 3(a) and the corresponding dip in Fig. 3(b) at β = 6.3 show the location of the 3:1 MMR, where instabilities occur over longer time-scales.

Figure 3.

Similar to Fig. 2 but for three (panels a and b) and five (panels c and d) Pluto-mass moons.

For five Pluto-mass moons (in Figs 3c and d), there are very narrow ranges (4.8 ≲ β ≲ 5.1) for which the system can survive up to 10⁷ orbits of the innermost moon. As shown earlier for a system of three moons, the 2:1 MMR delineates a lower boundary and puts a limit on stability at β ∼ 4.5, but the outermost moon’s eccentricity grows to nearly 1 at β ∼ 5.5, where the perturbations from the Sun are driving its eccentricity to such high values. From Table 2, the outermost moon begins beyond the outer stability limit (Rosario-Franco et al. 2020) at β ∼ 5.5 (yellow dashed line in Fig. 3d). In between these boundary conditions (MMR overlap and perturbation from the Sun), a system of five Pluto-mass moons could be stable if the apsidal precession rate of the moons is similar enough to prevent orbital overlap or the moons remain out of phase in an MMR to avoid collision (e.g. the 3:2 resonance between Neptune and Pluto).

3.3 Case study 3: Luna-mass moons

Table 2 shows that the number of Luna-mass moons between the Roche radius and the Hill radius is limited to 4 due to the requirement of |$\beta _{\rm min} = 2\sqrt{3}\approx 3.5$|⁠. Thus, we perform simulations considering only three and four Luna-mass moons following a similar procedure as in Sections 3.1 and 3.2. The maximum eccentricity and system lifetime for three moons are shown in Figs 4(a) and (b), while Figs 4(c) and (d) illustrate the same measures for four Luna-mass moons. Three moons can maintain stable orbits for a narrow range in spacing, 4 < β < 6. An additional moon (Figs 4c and d) shows that the MMRs encountered by the outermost moon from the Sun and the overlap of MMRs between moons destabilize the entire system.

Figure 4.

Similar to Fig. 2 but for three (panels a and b) and four (panels c and d) Luna-mass moons.

Quarles & Lissauer (2018) explored planet packing in α Centauri, where the secondary star significantly perturbs the planetary system, and they showed that three-planet systems can be long-lived with an appropriately chosen orbital spacing and initial eccentricity of the planets. The system architecture of planets orbiting one star of a stellar binary is similar to our system of moons orbiting an Earth-mass planet, where the forced eccentricity and boundaries of orbital stability limit the total number of satellites. Due to the similarity in structure, we arrive at similar conclusions.

3.4 Log-linear fit of unstable lifetimes

To compare the system lifetimes of different orbital architectures in planet packing, previous studies (Chambers et al. 1996; Smith & Lissauer 2009; Obertas et al. 2017; Quarles & Lissauer 2018) employ a log-linear fit to the first transition region to stability (β ∼ 4–8). In Figs 2–4, we mark these approximately linear trends using cyan lines, where we measure the slope b and y-intercept c. Quarles & Lissauer (2018) showed that a constant β-shift to account for |$\beta _{\rm min}=2\sqrt{3}$| is necessary to ensure a fair comparison between simulated systems. As a result, we fit the unstable data points in the transition region to a log-linear function of the form

$$\begin{eqnarray*} \log _{10}\, t = b^\prime \beta ^\prime + c^\prime , \end{eqnarray*}$$

(8)

where the prime (′) coordinates refer to fits made with a shift in the x-axis (i.e. |$\beta ^\prime = \beta -2\sqrt{3}$|⁠). Consequently, this shift in the x-axis minimizes the correlation between the slope (b′) and the y-intercept (c′). This fit function is similar to the ones used by Quarles & Lissauer (2018) and Lissauer & Gavino (2021) and has allowed us to make a consistent comparison of the coefficients from the previous work.

Table 3 shows a general trend in the coefficients of the log-linear function for all moon types, where fewer moons (n) correspond to a steeper slope (b′) and a longer system lifetime (c′) at β_min. The slope of the fit indicates how the system lifetime changes as a function of the orbital spacing parameter β. For example, a steeper slope conveys that an increase in β significantly extends the system lifetime.

Table 3.

Coefficients for the log-linear fits (cyan coloured lines in Figs 2–4) using log₁₀(t) = b′β′ + c′ (equation 8). The primed values (m′ and b′) constitute the shift in β by |$2\sqrt{3}$|⁠.

Mass	n	Slope	y-intercept
		(b′)	(c′)
Ceres	3	1.37	2.21
	4	1.32	1.83
	5	1.31	1.59
	6	1.25	1.56
	7	1.24	1.53
	8	1.24	1.51
Pluto	3	6.45	0.50
	4	5.35	0.24
	5	5.01	0.20
Luna	3	12.44	1.96
	4	6.4	1.88
Earth-massa	3	0.996	2.234
	5	0.742	2.084
Earth-massb	3	1.68	1.799

^aA system of up to five Earth-mass planets orbiting α Cen B (Quarles & Lissauer 2018).

^bA hypothetical system, where three Earth-like planets orbit a Sun-like star (Lissauer & Gavino 2021).

For Ceres-mass moons, the slope only decreases by |$\sim 10{{\ \rm per\ cent}}$| as more moons are added. The slopes for Ceres-mass moons are more similar to values determined through packed three-planet systems around a single star (Lissauer & Gavino 2021) rather than within a binary like α Centauri AB (Quarles & Lissauer 2018). This is likely due to the minimal forcing from the Sun, especially when the moons occupy a smaller portion of the host planet’s Hill radius. Table 2 suggests that more Ceres-mass moons could stably orbit the host planet, but the moons are perturbed and scattered by the neighbouring moons (or the star) due to their lower inertia. Hence, the system of Ceres-mass moons is largely unstable for β spacing ≤ 6.5. The decreasing slopes also indicate that the orbital spacing β must be increased, for increasing n, in order to avoid orbital crossings and maintain stability.

The other massive moons (Pluto-mass and Luna-mass) can absorb more internal (moon–moon) perturbation due to the increased inertia, or smaller changes to their angular momentum. Hence, the slopes for three-moon cases are much steeper compared to the Ceres-mass (6 times for Pluto-mass and 12 times for Luna-mass). In addition, pairs of more massive moons have a wider mutual Hill radius, which reduces the number of moons that can fit within the stability boundary. Therefore, the number of stable moons for Ceres-mass is ≤ 8, for Pluto-mass is ≤ 5, and for the Luna-mass is ≤ 4. The decreasing slopes for Pluto-mass and Luna-mass moons enforce our previous conclusion that the β must be increased, for increasing n, to avoid the orbital crossings.

3.5 Analysis of MMR using MEGNO maps

We explore the potential routes to instability in systems of multiple moons by using the chaos indicator MEGNO (mean exponential growth of nearby orbits 〈Y〉) maps. The MEGNO criterion is generally used to distinguish between chaotic, periodic, and aperiodic orbits within a phase space. Analysis of many-body systems using MEGNO was originally developed by Cincotta & Simó (1999, 2000) and Cincotta, Giordano & Simó (2003) to identify potential instabilities due to resonance overlap more efficiently in a short numerical integration period. It is capable of detecting high-order resonances (see e.g. Satyal et al. 2014 for a MEGNO map of a circumprimary planet displaying 39:2 MMR in a binary system) due to its sensitivity to unstable orbits and is a global indicator of dynamical changes in any Hamiltonian system. In our case, we are exploiting this particular nature of MEGNO, which can reveal fine resonance structures in a phase space to display the unstable and chaotic regions.

We have limited the phase space up to e_o = 0.3 because (for higher e_o) it is mostly chaotic region induced by the secular evolution of the eccentricities (see e.g. fig. 8 of Tamayo et al. 2021).

Figure 5.

MEGNO map for a system with three Ceres-mass moons simulated for ∼2 million orbits of the innermost moon. Panels (a) and (b) show the measure of MEGNO 〈Y〉 for the modified system (i.e. Earth and moons only) and full system (i.e. Sun–Earth–moons), respectively. MEGNO 〈Y〉 values equal to 2 indicate periodic orbits, while higher values indicate chaotic and potentially unstable orbits. Panel (c) shows the variation in β between the first two moons that is attained over a simulation. The green vertical lines represent the estimated locations of the MMRs (Murray & Dermott 2000). The MMR offset between the expected locations and observed numerical solutions shows the influence of the Sun on the satellite systems.

Figure 6.

Similar to Fig. 5(b), but for a system with three (a) Pluto-mass and (b) Luna-mass moons.

Figure 7.

Orbital evolution of three Ceres-mass moons with respect to their scaled distance from the host planet d/R_H (panels a and e), eccentricity e (panels b and f), orbital spacing β (panels c and g), and the 5:3 resonant angle ϕ between pairs of moons (panels d and h). The systems are simulated for 25 yr with an initial orbital spacing β = 7.26 (i.e. the resonance location for the 5:3 MMR in Fig. 6). The orbital evolution on the left-hand side includes the Sun in the N-body simulation, whereas the panels on the right-hand side exclude the Sun.

Figure 8.

Outward migration within a CTL tidal secular model of a Sun–Earth–moon system that varies the (a) time lag or (b) assumed satellite mass. In panel (a), the satellite mass is 3 m_Luna, while the time lag varies (colour-coded). In panel (b), the satellite mass is varied (three times Ceres-, Pluto-, and Luna-mass) and the time lag Δt is 100 s. The grey region marks the unstable region, past the stability limit, for prograde orbits.

Fig. 5(a) displays MEGNO values (colour-coded) calculated for the outermost moon in a system with an Earth-mass planet and three Ceres-mass moons (ignoring the Sun) with a wide variation in the initial β spacing and eccentricity (of the outermost moon) in the e_o–β phase space. MEGNO values 〈Y〉 = 2 correspond to periodic (and presumably stable) orbits, whereas 〈Y〉 = 6 indicates initial conditions that drive the system into chaos. Other quasi-periodic orbits can exist between these extremes (purple-orange). The outermost moon’s initial eccentricity is varied from 0 to 0.3 and the spacing β between all the moons ranges from 3.5 to 9.5. Most regions in the chaotic regions of the phase space (yellow) are unstable, where the outermost moon has a high probability of collision with the host planet or one of the moons. The existence of the chaotic regions is largely due to MMR overlap or secular evolution of the eccentricity (e.g. Wisdom 1980; Mudryk & Wu 2006; Quillen 2011; Laskar & Petit 2017; Hadden & Lithwick 2018; Petit et al. 2020; Tamayo et al. 2021). The stable orbits (black) begin when β ≳ 6, and for lower initial eccentricity. For one initial condition (e_o = 0), the map compliments the analysis presented in Section 3.1, Fig. 2(a), where the orbits are stable for the full integration period when a higher β spacing is selected. Figs 5(b) and (c) both explore a similar phase space, but now including the Sun in the simulations.

Multiple MMRs are V-shaped (because their libration width increases with increasing eccentricity) in the map, where the 8:5, 5:3, 7:4, and 9:5 are some of the prominent MMRs between moons and are marked (dashed green lines). Fig. 5(b) demonstrates how the solar perturbations affect the initial conditions that are affected by the moon–moon MMRs. As a result, there is an offset observed in the MEGNO phase space when comparing Figs 5(a) and (b). For example, Fig. 5(a) shows the 5:3 MMR located at β = 7.2, but the MMR structure is located at β = 7.4 in Fig. 5(b). This shift is accounted for by the secular perturbations on the moons by the Sun, which causes the orbital spacing of the moons to change over time.

A moon can then evolve into and out of a nearby MMR. Fig. 5(c) illustrates the shifts in orbital spacing between the two inner moons through Δβ₁₂, which signifies the maximum change in the orbital spacing. In and around the observed MMRs, the shift in the Δβ is about 0.1, which accounts for the differences in the location of MMRs in Figs 5(a) and (b). Also observed in Fig. 5(c) is that the Δβ values higher than 0.3 indicate unstable orbits and those values that are less than 0.3 are stable. It is also apparent from the map that the orbits with Δβ less than 0.2 are not heavily influenced by any MMRs.

A similar phase space is explored for larger moons, but focusing on three-moon systems with both the host planet and the Sun included in Figs 6(a) (Pluto-mass moons) and (b) (Luna-mass moons). Fig. 6(a) shows that the outermost Pluto-mass moon can start with relatively large eccentricity as long as the orbital spacing is large enough (β ≳ 6.5), while an equivalent system with Luna-mass moons is more limited in terms of its initial orbital spacing (Fig. 6b). Both systems require larger exchanges of angular momentum to perturb their orbits, as compared to the Ceres-moon case (Fig. 5b), which explains why a higher initial eccentricity can allow for stable, periodic orbits (black regions). Pluto-mass moons appear the most optimal as nearly half of the parameter space allows for periodic orbits. However, the MMR at β ∼ 6.4 may induce instabilities as a primordial moon system evolves outwards due to tidal interactions with the host planet. Also, at this resonance, the moons’ e_max are observed to evolve towards 1 (Fig. 3a).

3.6 Shifting the MMRs

The inclusion (or exclusion) of the Sun and its effect on the dynamics of the moons are viewed in global phase space maps (e_o versus β) using the MEGNO criterion (Fig. 5). For the case with three Ceres-mass moons, the 5:3 MMR is clearly observed at β = 7.26 (without the Sun, Fig. 5a) and slightly shifted higher at β = 7.4 (with the Sun, Fig. 5b). To visualize the time-series data of individual orbital elements, we use initial conditions (β = 7.26 and e_o = 0.01) from Fig. 5 and simulate the system for 25 yr (∼3 × 10⁶ orbits of the innermost moon). Fig. 7 shows the time-series evolution of the normalized orbital distance d/R_H, eccentricity e, orbital spacing β, and resonant angle ϕ_{5: 3} for each moon in both systems (with and without including the Sun). Note that the two (inner and middle) moons begin on circular orbits, while the outer moon begins with the eccentricity e_o = 0.01, to maintain consistency with Fig. 5. The normalized distance (d/R_H) is used instead of semimajor axis to better illustrate the correlated changes in distance with eccentricity as it affects the instantaneous measure of the orbital spacing.

In Fig. 5(b) (β = 7.26 and e_o = 0.01 with the Sun included), the value of MEGNO suggests that the orbits are periodic. This is confirmed in the evolution of each moon’s normalized distance (Fig. 7a) and eccentricity (Fig. 7b), where the colours refer to the inner (red), middle (blue), and outer (green) moons. The gravitational perturbation of the Sun forces a high-frequency variation in the relative distance of each moon, which underlies the variation in the orbital spacing β (Fig. 7c). This forcing also prevents the 5:3 MMR resonant angle of the inner pair (ϕ₁₂) or the outer pair (ϕ₂₃) of moons from librating (in Fig. 7d), even though the respective orbital spacing of each pair (β₁₂ in orange or β₂₃ in cyan) crosses the expected MMR location β = 7.26. In contrast, Figs 7(e)–(h) display a similar simulation, where the Sun is removed (i.e. ignoring its secular perturbation). In this case, the normalized distance and the eccentricity of each moon are chaotic, which is indicated by the corresponding MEGNO value from Fig. 5(b). The variation of each moon’s eccentricity (in Fig. 7f) can be two to four times larger as compared to Fig. 7(b) due to the eccentricity excitation from the 5:3 MMR. The orbital spacing of the inner and outer moon pairs varies on much slower time-scales in Fig. 7(g), which allows for moon pairs to evolve together into the 5:3 MMR with an average β ∼ 7.26. The switching between libration and circulation in Fig. 7(h) confirms the chaotic nature of this initial condition and the connection with the 5:3 MMR.

For an Earth-mass planet with a semimajor axis of 1 au, the Sun clearly contributes to the dynamics of the planet and its moons. The inclusion of the Sun in a system is necessary for our analysis. However, the chaotic dynamics that remains by removing the Sun is still applicable, where the relevant initial orbital spacing β must be shifted by ∼0.15 to β ∼ 7.4. Alternatively, moon systems of long-period Earth-mass planets could undergo similar chaotic variations (e.g. Pluto–Charon and its four moons).

3.7 Mass distribution and formation plausibility

Based on the number of moons that stably orbit the planet, we calculate the total mass distribution within the stability boundary for all three systems. The total mass of three Luna-mass, five Pluto-mass, and eight Ceres-mass moons is equivalent to 1.11 × 10⁻⁷, 3.30 × 10⁻⁸, and 3.60 × 10⁻⁹ M_⊙, respectively. Despite the small differences in the semimajor axis of the innermost moon a₁, the order-of-magnitude differences in mass (see Table 1) are preserved. Since the mass is distributed roughly in the same surface area between the inner (i.e. Roche limit) and outer stability boundaries, the surface mass density for Ceres (σ_Ceres) is approximately 7.2 × 10⁻⁵ M_⊙ au⁻², equivalent to 637 g cm⁻². The surface mass densities for Pluto and Luna are 10σ_Ceres and 30σ_Ceres, respectively.

Our work concentrates on the maximum number of satellites, with different masses, that can stably orbit around an Earth-like planet. Whether more than one satellite can form around an Earth-like planet is beyond the scope of this work. However, we can provide some assessment on the plausibility of formation based upon our results of the total mass surface density for the maximum number of moons. Moons can form around a giant planet (e.g. Jupiter and Saturn) from a gaseous circumplanetary disc during the last stages of planet formation as has been shown for the Galilean moons (Pollack & Bodenheimer 1989; Canup & Ward 2006). But, there is not a known minimum planetary mass threshold at which a circumplanetary disc can form and evolve into moons (Ayliffe & Bate 2009). Large moons are expected to arise, primarily due to giant impacts, where Nakajima et al. (2022) suggests that impact-induced large moons are more likely to form around rocky planets whose radius is smaller than 1.6 R_⊕. Moreover, numerical models using smoothed particle hydrodynamics (SPH) have shown that the mass surface density of moon-forming discs can reach ∼10⁷ g cm⁻² only at a few Earth radii and substantially spread over the lifetime of the disc (Nakajima & Stevenson 2014; Nakajima et al. 2022). Therefore, it is at least plausible that multiple moons could form around Earth-like planets. Further study or confirmation of the current exomoon candidates (e.g. Kepler 1625b-i, Teachey & Kipping 2018; Kepler 1708b-i, Kipping et al. 2022) could shed more light on these hypotheses.

3.8 Effects of tides

The outward migration of satellites through tidal interactions modifies their potential lifetime (Barnes & O’Brien 2002; Sucerquia et al. 2019; Lainey et al. 2020). In the Solar system, Charnoz, Salmon & Crida (2010) suggested that the population of small moons that orbit just outside Saturn’s rings could have originated at the edge of the main rings and tidally migrated outwards. To obtain a complete picture of the orbital stability of exomoons, it is necessary to consider the contribution of planetary and stellar tides. We apply a secular constant time lag (CTL) tidal model (Hut 1981; Leconte et al. 2010) and evaluate the migration time-scales of moons assuming that moon formation readily occurs near the host planet’s Roche limit. The lifetime of a moon system can be reduced as the outermost moon migrates towards the outer stability limit (i.e. ∼0.40 R_H; Rosario-Franco et al. 2020), where this will depend on the mass of the satellite (or moon–planet mass ratio) and the assumed time lag Δt for the tidal dissipation. The secular model calculates the changes to the orbital elements of both the host planet and its moon through the respective semimajor axes (a_p and a_sat), eccentricities (e_p and e_sat), and mean motion (n_p and n_sat) averaged over an orbit. The model is scaled by the tidal Love number k₂ and the time lag Δt, where the latter is proportional to (nQ)⁻¹ in the constant phase lag (i.e. constant Q) tidal models (Leconte et al. 2010; Piro 2018).

We consider two scenarios for three-moon systems: (a) keep the satellite mass fixed at 3 m_Luna, while evaluating a range of CTL values, and (b) keep the time lag fixed at Δt = 100 s, while evaluating a range in satellite mass (Ceres-, Pluto-, and Luna-mass). The CTL Δt is varied between simulations from 10 to 600 s on a logarithmic scale. The moon–moon interactions are ignored since we are applying a secular model, where the total mass of the moon system is combined into the innermost moon. This represents a conservative estimate because the innermost moon would be migrating outwards more slowly in reality as it would be 1/3 of our prescribed mass. The Earth-mass host planet is assumed to begin with a rotation period of 5 h, which is consistent with expectations from terrestrial planet formation (Kokubo & Ida 2007), and we are interested in a 10¹⁰ yr time-scale (i.e. the main-sequence lifetime of a G dwarf).

In our first scenario, shown in Fig. 8(a), the mass of the satellite is 3 m_Luna and the CTL is varied over a range that corresponds to a very low dissipation (10 s; red) up to a very high dissipation (600 s; lavender). In all cases, the satellite’s outward migration stalls at ∼0.1 R_H as the moon’s orbital period synchronizes with the planet’s rotation period. Assuming that the moons migrate outwards together maintaining an orbital spacing of β = 4, then the outermost moon would migrate beyond the outer stability limit (i.e. a₃ ∼ 0.55 R_H) and thereby reduce the number of stable moons by one.

The effect of tidal migration is more dire for lower mass moons because they can migrate closer to the outer stability limit (Fig. 8b) and require a larger orbital spacing (see Figs 5b and 6a). This combination of circumstances will likely cause at least one moon to scatter and/or migrate past the outer stability limit. Therefore, outward tidal migration will likely reduce the number of moons orbiting an Earth-mass planet by at least one in three-moon systems and likely more within moon systems of higher multiplicity. Further consideration of tides is beyond the scope of our work, where others could explore the effects of differential tides on outward migration that is similar to models for Saturn’s moons (Crida & Charnoz 2012; Ćuk, Dones & Nesvorný 2016).

4 CONCLUSIONS

Through N-body simulations, we investigate the potential for systems of three to nine moons orbiting an Earth-mass planet and a solar-mass star. The moons vary in mass, but are analogous to Ceres, Pluto, and our Moon (Luna). Systems of multiple moons are inherently constrained by the inner Roche limit and the outer stability limit, which can be also scaled by a planet’s Hill radius. Scaling by the Hill radius allows our work to be generalized beyond an exact Earth–Sun analogue for the primary bodies, because the Hill radius incorporates potential changes in the planetary semimajor axis, eccentricity, and mass, in addition to the stellar mass. Each moon system begins on a circular and coplanar orbit, where the initial orbital phase is selected through the golden ratio following planet-packing studies (Smith & Lissauer 2009; Quarles & Lissauer 2018; Lissauer & Gavino 2021). We find using N-body simulations that 7 ± 1 Ceres-mass moons could stably orbit an Earth-mass planet at 1 au from a Sun-like star. If the moons are more massive (Pluto- or Luna-mass), then the number of moons with stable orbits reduces to 4 ± 1 and 3 ± 1, respectively. Outward tidal migration will likely modify these estimates by at least one moon, where additional moons could be lost through scattering, collisions, or simply migrating beyond the outer stability limit.

The orbital spacing between each moon is measured using a dimensionless parameter β, which is the distance between two neighbouring moons divided by their mutual Hill radius. The maximum number of moons produces a minimum in the orbital spacing, where we find β = 6, 4.5, and 3.5 for Ceres-, Pluto-, and Luna-mass moons, respectively. The potential stability for these moon systems depends on their proximity to MMRs between adjacent pairs of moons. The location of the MMRs is estimated using the chaos indicator MEGNO, which shows a shift of ∼0.15 in β from the expected location due to perturbations on the moon system from the Sun. We show that a moon can behave chaotically (i.e. periodic switching between circulation and libration in the 5:3 resonant angle) when starting within the libration zone of the 5:3 MMR and ignoring the gravitational solar perturbations.

Planet-packing studies have used a best-fitting slope from a log-linear model to compare the changes in stability due to planet multiplicity (Chambers et al. 1996; Smith & Lissauer 2009; Pu & Wu 2015; Obertas et al. 2017; Quarles & Lissauer 2018). We employ a similar technique, but in shifted coordinates so that the y-intercept occurs at |$\beta = 2\sqrt{3}$| (e.g. Quarles & Lissauer 2018; Lissauer & Gavino 2021). From these measurements, we find that Ceres-mass moons have a slope in the range 1.24–1.37, which is inversely correlated with the number of moons. These slopes are less than the expected values for systems of three Earth-mass planets orbiting a solar-mass star and greater than the more extreme case where three planets are orbiting α Centauri B, as stellar binary. Pluto- and Luna-mass moons have a much steeper slope because they have a larger mutual Hill radius, which drastically limits the potential orbital spacing between moons (see Table 2).

We use the mass surface density within the stability boundary limits to determine whether systems of multiple moons are at least plausible. Nakajima et al. (2022) showed that the surface density of a moon-forming disc can reach ∼10⁷ g cm⁻², which is much higher than the mass surface density of three Luna-mass moons (∼10⁴ g cm⁻²). The mass surface density for five Pluto-mass moons is smaller by a factor of 3, while it decreases by a factor of 30 for eight Ceres-mass moons. It appears plausible that multiple moons could form, but further study using SPH simulations would be necessary and is beyond the scope of our work.

Our N-body simulations are mostly limited to only 10⁷ orbits of the innermost moon, which is ∼3 × 10⁴ yr. The long-term evolution of moon systems will be determined by the outward migration of the moons due to tidal interactions with the host planet. We evaluate this possibility using a CTL secular model (Hut 1981; Leconte et al. 2010; Quarles et al. 2021), where outward tidal migration becomes significant at long time-scales (∼10⁸–10¹⁰). As a result, the number of moons that can stably orbit an Earth-mass planet is reduced by one. In the case of Luna-mass moons, only one moon is lost, where systems of less massive moons can have more significant losses due to the relative ease of scattering events and the final migration distance of the innermost moon.

Detecting multiple-moon systems orbiting other stars is currently out of reach, where there are only a couple of exomoon candidates using the photometric detection method (Teachey & Kipping 2018; Kipping et al. 2022). Recent observations from Atacama Large Millimeter/submillimeter Array (Benisty et al. 2021) suggest that a moon-forming disc exists around PDS 70c, which points to the potential for long-wave observations or direct imaging (Vanderburg, Rappaport & Mayo 2018). From this work, the dynamical stability of moon systems limits the existence of exomoons for Earth analogues in their respective habitable zones, where confirmation from future observations is needed.

ACKNOWLEDGEMENTS

The authors appreciate the constructive comments and feedback from the referee. MRF acknowledges support from the NRAO Gröte Reber Fellowship and the Louis Stokes Alliance for Minority Participation Bridge Program at the University of Texas at Arlington. This research was supported in part through research cyberinfrastructure resources and services provided by the Partnership for an Advanced Computing Environment (PACE) at the Georgia Institute of Technology.

DATA AVAILABILITY

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